Method for performing parallel magnetic resonance imaging

ABSTRACT

A method of parallel magnetic resonance imaging of a body, comprising:—acquiring a set of elementary magnetic resonance images of said body from respective receiving antennas having known or estimated sensibility maps and noise covariance matrices, said elementary images being under-sampled in k-space; and performing regularized reconstruction of a magnetic resonance image of said body; wherein said step of performing regularized reconstruction of a magnetic resonance image is unsupervised and carried out in a discrete frame space. A method of performing dynamical and parallel magnetic resonance imaging of a body, comprising:—acquiring a set of time series of elementary magnetic resonance images of said body from respective receiving antennas having known or estimated sensibility maps and noise covariance matrices, said elementary images being under-sampled in k-space; and performing regularized reconstruction of a time series of magnetic resonance images of said body.

FIELD

The invention relates to a method for performing parallel magneticresonance imaging (pMRI) of a body, including parallel, dynamical(time-resolved) magnetic resonance imaging, such as functional MRI(fMRI).

BACKGROUND

Reducing global acquisition time is of main interest in medical magneticresonance imaging (MRI), and even more when dynamic imaging such as fMRIis concerned. Actually, a short acquisition time allows improving thespatial/temporal resolution of acquired fMRI data, which leads to a moreefficient statistical analysis. In addition, by reducing the globalimaging time, some additional artifacts caused by the patient motion canbe avoided. For this reason, parallel imaging systems have beendeveloped: multiple receiver surface coils with complementarysensitivity profiles located around the underlying object or body areemployed to simultaneously collect in the frequency domain (i.e. theso-called k-space), data sampled at a rate R times lower than theNyquist sampling rate along at least one spatial direction, i.e. thephase encoding one; R is usually called the “reduction factor”.Therefore, the total acquisition time is R times shorter than withconventional non parallel imaging. A reconstruction step is thenperformed to build a full Field of View (FOV) image by unfolding theundersampled “elementary” images acquired by the individual receivers.This reconstruction is a challenging task because of the low Signal toNoise Ratio (SNR) in parallel MRI (pMRI) caused by aliasing artifactsrelated to the undersampling rate, those caused by noise during theacquisition process and also the presence of errors in the estimation ofcoil sensitivity maps.

The Simultaneous Acquisition of Spatial Harmonics (SMASH) [Sodickson etal., 1997] was the first reconstruction method, operating in the k-spacedomain. It uses a linear combination of pre-estimated coil sensitivitymaps to generate the missing phase encoding steps.

Some other k-space based reconstruction techniques have also beenproposed like GRAPPA (Generalized Autocalibrating Partially ParallelAcquisitions) [Griswold et al., 2002], and SENSE (Sensitivity Encoding)[Pruessmann et al., 1999]. SENSE is a two-step procedure relying firston a reconstruction of reduced FOV images and second on a spatialunfolding technique, which amounts to a weighted least squaresestimation. This technique requires a precise estimation of coilsensitivity maps using a reference scan (usually a 2D Gradient-Echo(GRE)). It is presently the most frequently employed pMRI technique,applied in particular to brain and cardiac imaging.

For a general overview of reconstruction methods in pMRI see [Hoge etal., 2005].

SENSE is often supposed to achieve an exact reconstruction in the caseof noiseless data and perfect coil sensitivity maps knowledge, which isalso true for all above mentioned methods. However, in practice, thepresence of noise in the data and inaccuracies in the estimation of coilsensitivity maps are unavoidable and make the reconstruction problemill-conditioned.

As image reconstruction is an ill-posed inverse problem, regularizationtechniques are commonly applied to better estimate the full FOV image.Most of these techniques operate in the image domain; in particular,this is the case for Tikhonov regularization [Ying et al., 2004], whichuses a quadratic penalty term either to promote smoothness constraintsor to account for the squared difference between the reconstructed imageand an a priori reference image. Despite the use of regularization,however, high reduction factors (exceeding a value of R=2) are generallyconsidered unfeasible when low magnetic field intensities (up to 1.5Tesla) are used, since the reconstructed images are affected by severealiasing artifacts.

In [Chaâri et al. 2008] and [Chaâri et al. 2009] the present inventorshave described a method of performing regularized image reconstructionin parallel MRI, using a wavelet-based regularization scheme, allowingto increase the reduction factor R.

The present invention aims at providing several improvements of saidmethod, including extending it to dynamical imaging (e.g. fMRI) andmaking it fully or partially auto-calibrated (or “unsupervised”).

SUMMARY

An object of the present invention is then a method of parallel magneticresonance imaging of a body, comprising:

-   -   acquiring a set of elementary magnetic resonance images of said        body from respective receiving antennas having known or        estimated sensibility maps and noise covariance matrices, said        elementary images being under-sampled in k-space; and    -   performing regularized reconstruction of a magnetic resonance        image of said body;

wherein said step of performing regularized reconstruction of a magneticresonance image is carried out in a discrete frame space by minimizing acost function comprising:

-   -   an error term, representative of a likelihood of a reconstructed        image, given said acquired elementary images; and    -   a frame penalty term, representative of a deviation between an        actual statistical distribution of frame coefficients of said        reconstructed image and an a priori distribution of said        coefficients;

said a priori distributions of the frame coefficients of thereconstructed image being estimated on the basis of an auxiliarymagnetic resonance image of said body.

According to different embodiments of said method:

-   -   Said error term can be representative of a neg-log-likelihood of        said reconstructed image, given said acquired elementary images.    -   Said step of performing regularized reconstruction of a magnetic        resonance image of said body can be carried out by maximizing,        in said frame space, a full posterior distribution of a set of        frame coefficients defining an image of the body, given said        acquired elementary magnetic resonance images and said a priori        distribution of the frame coefficients.    -   Said auxiliary magnetic resonance image of said body can be        reconstructed from said acquired elementary magnetic resonance        images. More particularly, said auxiliary magnetic resonance        image of said body can be reconstructed from said acquired        elementary magnetic resonance images using a SENSitivity        Encoding—SENSE—algorithm. Even more particularly, said auxiliary        magnetic resonance image of said body can be reconstructed from        said acquired elementary magnetic resonance images using an        algorithm chosen between: an unregularized SENSE algorithm; a        SENSE algorithm regularized in image space; and a SENSE        algorithm regularized in k-space.    -   A generalized Gauss-Laplace a priori statistical distribution of        said frame coefficients can assumed, and parameters of said        distribution can be estimated on the basis of said auxiliary        magnetic resonance image of said body, using a        maximum-likelihood or a posterior mean estimator.    -   Said error term can be a quadratic mean error term.    -   Said acquired elementary images can be three-dimensional images,        and said step of performing regularized reconstruction of a        magnetic resonance image can be carried out in a discrete        three-dimensional frame space. More particularly, said acquired        three-dimensional elementary images can be obtained by stacking        bi-dimensional elementary images of slices of the object to be        imaged.    -   Said step of performing regularized reconstruction of a magnetic        resonance image can be based on a redundant wavelet frame        representation. Alternatively, said step of performing        regularized reconstruction of a magnetic resonance can be based        on a non-redundant wavelet representation.    -   Said cost function can also comprise at least one spatial domain        penalty term chosen between: a total variation norm of the        reconstructed image; and a convex constraint.

Another object of the present invention is a method of performingdynamical and parallel magnetic resonance imaging of a body, comprising:

-   -   acquiring a set of time series of elementary magnetic resonance        images of said body from respective receiving antennas having        known or estimated sensibility maps and noise covariance        matrices, said elementary images being under-sampled in k-space;        and    -   performing regularized reconstruction of a time series of        magnetic resonance images of said body;

wherein said step of performing regularized reconstruction of a timeseries of elementary magnetic resonance images is carried out byminimizing a cost function comprising:

-   -   an error term, representative of a likelihood of each        reconstructed image, given corresponding acquired elementary        images; and    -   a temporal penalty term, representative of a pixel-by-pixel or        voxel-by-voxel difference between consecutive image of the        series.

According to different embodiments of said method:

-   -   Said temporal penalty term can be based on an edge-preserving        function, more particularly a convex edge-preserving function        and even more particularly on an L_(p) norm with p. 1 and        preferably 1≤p<1.5    -   Said temporal penalty term can be given by the sum of a first        partial temporal penalty term and a second partial temporal        penalty term, wherein: the first partial temporal penalty term        is representative of pixel-by-pixel or voxel-by-voxel        differences between each even-numbered image of the series and a        preceding odd-numbered image; and the second partial temporal        penalty term is representative of pixel-by-pixel or        voxel-by-voxel differences between each odd-numbered image of        the series and a preceding even-numbered image; said cost        function being minimized by using proximity operators for said        first and second partial temporal penalty terms.    -   Said step of performing regularized reconstruction of a magnetic        resonance image can be carried out in a discrete frame space,        and said cost function can also comprise a frame penalty term,        representative of a deviation between statistical distributions        of frame coefficients of each reconstructed image and an a        priori distributions of said coefficients; said a priori        distributions of the frame coefficients of the reconstructed        images being estimated on the basis of an auxiliary magnetic        resonance image of said body. In this case, said elementary        images can be three-dimensional images, and said discrete frame        space can be a discrete three-dimensional frame space.    -   The method can further comprise a step of automatically        determining, using a maximum-likelihood estimator; a weighting        parameter of said temporal penalty term. More precisely, the        method can comprise estimating said weighting parameter of the        temporal penalty term for each pixel or voxel, or set of        neighboring pixels or voxels, of the image to be reconstructed.        Moreover, said temporal penalty term can be based on an L_(p)        norm and said weighting parameter of the temporal penalty term        and the parameter p can be jointly determined using said        maximum-likelihood estimator. In particular, said weighting        parameter of the temporal penalty term and the parameter p can        be jointly determined using said maximum-likelihood estimator        under the constraint p≥1.    -   Said error term can depend on geometrical parameters defining a        rigid transformation of each of said elementary magnetic        resonance images with respect to an elementary magnetic        resonance image taken as a reference, and wherein said step of        performing regularized reconstruction of a time series of        elementary magnetic resonance images is carried out by        minimizing said function also with respect to said geometrical        parameters.    -   Said elementary images can be acquired by echoplanar imaging.    -   Said elementary images can be under-sampled with a reduction        factor higher or equal to 4.

BRIEF DESCRIPTION OF THE DRAWINGS

Additional features and advantages of the present invention will becomeapparent from the subsequent description, taken in conjunction with theaccompanying drawings, which show:

FIG. 1, empirical histograms of real and imaginary parts of waveletcoefficients of a magnetic resonance image of a human brain;

FIGS. 2 and 3, nine anatomical axial slices of a human brain obtainedusing the Tikhonov-regularized SENSE pMRI reconstruction method knownfrom prior art, with R=2 and R=4 respectively;

FIG. 4, an anatomical axial slice of the same human brain obtained usingthe TV-regularized SENSE pMRI reconstruction method known from priorart, with R=2 (left panel) and R=4 (right panel);

FIGS. 5 and 6, nine anatomical axial slices of the same human brainobtained using an autocalibrated 2D wavelet transform-basedregularization scheme (UWR-SENSE) according to an embodiment of theinvention, with R=2 and R=4 respectively;

FIGS. 7 and 8, nine anatomical axial slices of the same human brainobtained using a constrained, autocalibrated 2D wavelet transform-basedregularization scheme (CWR-SENSE) according to another embodiment of theinvention, with R=2 and R=4 respectively;

FIG. 9, nine anatomical axial slices of the same human brain obtainedusing an autocalibrated combined wavelet-total variation regularizationscheme according to another embodiment of the invention, with R=4 andusing image decomposition on a non-redundant orthonormal wavelet basis;

FIG. 10, nine anatomical axial slices of the same human brain obtainedusing an autocalibrated combined wavelet-total variation regularizationscheme according to another embodiment of the invention, with R=4 andusing image decomposition on a redundant wavelet frame constituted bythe union of two orthonormal bases;

FIG. 11, three anatomical axial slices of the same human brain obtainedusing: conventional non-parellel MRI (top row), mSENSE parallel MRIreconstruction with R=4 (middle row) and an autocalibrated 3D wavelettransform-based regularization scheme (3D-UWR-SENSE) according toanother embodiment of the invention (bottom row);

FIG. 12, three anatomical sagittal slices of the same human brainobtained using: conventional non-parellel MRI (top row), mSENSE parallelMRI reconstruction with R=4 (middle row) and an autocalibrated 3Dwavelet transform-based regularization scheme (3D-UWR-SENSE) accordingto another embodiment of the invention (bottom row);

FIG. 13, three anatomical axial slices of the same human brain (upperrow) and magnified details thereof (bottom row) obtained using: TVregularization (left); combined wavelet-total variation regularizationscheme according to an embodiment of the invention using imagedecomposition on a non-redundant wavelet basis (center) and combinedwavelet-total variation regularization scheme according to an embodimentof the invention using image decomposition on a redundant wavelet frameconstituted by the union of two orthonormal bases (right);

FIG. 14, Axial, Coronal and Sagittal slices of a human brain, acquiredusing Ecoplanar (EPI) fMRI and reconstructed using the mSENSE method,known from prior art, and an autocalibrated 2D wavelet transform-basedregularization scheme (4D-UWR-SENSE) according to an embodiment of theinvention;

FIG. 15, subject-level student-t maps, superimposed to anatomicalimaging, of the aC-aS contrast detected using EPI fMRI; data have beenreconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE methods,respectively, with R=2 (top of the figure) and R=4 (bottom of thefigure); sagittal, coronal and axial views are displayed;

FIG. 16, subject-level student-t maps of the aC-aS contrast for twosubjects, reconstructed at R=2 using the mSENSE, UWR-SENSE and4D-UWR-SENSE methods, respectively; sagittal, coronal and axial viewsare displayed;

FIG. 17, subject-level student-t maps, superimposed to anatomicalimaging, of the Lc-Rc contrast detected using EPI fMRI; data have beenreconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE methods,respectively, with R=2 (top of the figure) and R=4 (bottom of thefigure); sagittal, coronal and axial views are displayed;

FIG. 18, subject-level student-t maps of the Lc-Rc contrast for twosubjects, reconstructed at R=4 using the mSENSE, UWR-SENSE and4D-UWR-SENSE methods, respectively; sagittal, coronal and axial viewsare displayed;

FIG. 19, group-level student-t maps for the aC-aS contrast, where datahave been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSE forR=2 and R=4; sagittal, coronal and axial views are displayed; and

FIG. 20, group-level student-t maps for the Lc-Rc contrast, where datahave been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSE forR=2 and R=4; sagittal, coronal and axial views are displayed.

DETAILED DESCRIPTION

Before describing the invention in detail it will be necessary to recallsome basic facts about pMRI (and, in particular, SENSE).

MRI is an imaging technique that can proceed either in 2D(two-dimensional) or directly in 3D (three-dimensional) depending on theinvolved RF pulse design. In the 2D case, a volume is covered usingadjacent slices. PMRI may adapt to both situations since this methodmakes the k-space scanning faster, whatever its definition (i.e. in 2 or3D). For the sake of simplicity, we only present the 2D case in thefollowing.

An array of C coils is employed to measure the spin density p into theobject under investigation (e.g. a brain, in the head of a patient).Image acquisition is based on specific imaging sequences; in exemplaryembodiments of the invention, anatomical MRI is based on the 3D MPRAGEsequence while functional MRI involves is performed using 2D echoplanarimaging (EPI). The following description will be focused on the 2D case;then, the signal d_(c) received by each coil l (1≤c≤C) is the Fouriertransform of the desired 2D field p weighted by the coil sensitivityprofile s_(c), evaluated at some locations k=(k_(y), k_(x))^(T) in thek-space (the apex ^(T) means transposition). The received signal {tildeover (d)}_(c) is therefore defined by the sampling scheme:{tilde over (d)} _(c)(k _(r))=∫ρ(r)s _(c)(r)e ^(−ι2πk) ^(r) ^(T) ^(r)dr+ñ _(c)(k _(r))  (1)where ñ_(c) is a realization of an Additive White Gaussian Noise (AWGN)and r=(y, x)^(T) is the spatial position in the image domain. For thesake of simplicity, a Cartesian coordinate system is generally adopted.In its simplest form, SENSE imaging amounts to solving a one-dimensionalinversion problem due to the separability of the Fourier transform. LetΔy=Y/R be the sampling period where Y is the size of the field of view(FOV) along the phase encoding direction, let y be the position in theimage domain along the same direction, x the position in the imagedomain along the frequency encoding direction and R≤L the reductionfactor. A 2D inverse Fourier transform allows recovering the measuredsignal in the spatial domain. By accounting for the undersampling of thek-space by R, (1) can be re-expressed in the following matrix form:d(r)=S(r)ρ(r)+n(r)  (2)where:

$\begin{matrix}{{{{S(r)}\overset{\Delta}{=}{\begin{bmatrix}{s_{1}\left( {x,y} \right)} & \ldots & {s_{1}\left( {x,{y + {\left( {R - 1} \right)\Delta\; y}}} \right)} \\\vdots & \vdots & \vdots \\{s_{C}\left( {x,y} \right)} & \ldots & {s_{C}\left( {x,{y + {\left( {R - 1} \right)\Delta\; y}}} \right)}\end{bmatrix} \in {\mathbb{C}}^{C \times R}}}{d(r)}\overset{\Delta}{=}{{\begin{bmatrix}{d_{1}\left( {x,y} \right)} \\{d_{2}\left( {x,y} \right)} \\\vdots \\{d_{C}\left( {x,y} \right)}\end{bmatrix} \in {{\mathbb{C}}^{C}.{\overset{\_}{\rho}(r)}}}\overset{\Delta}{=}{\begin{bmatrix}{\overset{\_}{\rho}\left( {x,y} \right)} \\{\overset{\_}{\rho}\left( {x,{y + {\Delta\; y}}} \right)} \\\vdots \\{\overset{\_}{\rho}\left( {x,{y + {\left( {R - 1} \right)\Delta\; y}}} \right)}\end{bmatrix} \in {\mathbb{C}}^{R}}}}{and}{{n(r)}\overset{\Delta}{=}{{\begin{bmatrix}{n_{1}\left( {x,y} \right)} \\{n_{2}\left( {x,y} \right)} \\\vdots \\{n_{C}\left( {x,y} \right)}\end{bmatrix} \in {\mathbb{C}}^{C}}..}}} & (3)\end{matrix}$

In equation (2), (n(r))_(r) is a sequence of circular zero-mean Gaussiancomplex-valued random vectors. These noise vectors are i.i.d.(independent, identically distributed) and spatially independent withcovariance matrix Ψ of size C×C. In practice, Ψ is estimated byacquiring data from all coils without radio frequency pulses, and itsgeneric entry Ψ(c₁,c₂) corresponding to the covariance between the twocoils c₁ and c₂ is given by:

$\begin{matrix}{{{\Psi\left( {c_{1},c_{2}} \right)} = {\frac{1}{Y \times X}{\sum\limits_{({y,x})}{{{\underset{\_}{d}}_{c_{1}}\left( {y,x} \right)}{{\underset{\_}{d}}_{c_{2}}^{*}\left( {y,x} \right)}}}}},{\forall{\left( {c_{1},c_{2}} \right) \in {\left\{ {1,\ldots\mspace{14mu},C} \right\}^{2}.}}}} & (4)\end{matrix}$where (□)* stands for the complex conjugate. Note that a statisticalmodel can be assumed for matrix such as the statistical independencebetween coils or a nearest-neighbors statistical dependence model. Inthe first case, matrix Ψ becomes diagonal and the coil-dependentvariances σ_(n) _(c) ² can be estimated using the Mean AbsoluteDeviation (MAD) technique described by [Donoho, 1995], which relies onan additive Gaussian white noise (AGWN) assumption. In the simplestcase, Ψ can be simply be taken equal to the identity matrix. In thepresence of non-zero off-diagonal terms, the empirical covariance givenin Eq (4) is used to estimate them.

Therefore, the reconstruction step consists in inverting (2) andrecovering from d(r) at spatial positions r=(y, x)^(T). Note that thedata (d_(c))_(1≤|≤C) and the unknown image are complex-valued, although∥ is only considered for visualization purposes.

A simple reconstruction method also called the SENSE approach[Pruessmann et al., 1999], is based on the minimization of the WeightedLeast Squares (WLS) criterion.

The objective is to find a vector at each spatial location r such that:

$\begin{matrix}{{{\hat{\rho}}_{WLS}(r)} = {\arg\;{\underset{\rho{(r)}}{\;\min}{\mathcal{J}_{WLS}\left( {\rho(r)} \right)}}}} & (5)\end{matrix}$where

WLS(ρ(r))=∥d(r)−S(r)ρ(r)∥_(Ψ) ⁻¹ ², (·)^(H) stands for the transposedcomplex conjugate, (·)^(#) stands for the pseudo-inverse and defines anorm on.

As discussed above, this inverse problem is generally ill-posed andrequires regularization, e.g. Tikhonov regularization. Theregularization process typically consists in computing {circumflex over(ρ)}_(PWLS)(r) as the minimizer of the following Penalized WeightedLeast Squares (PWLS) criterion:

_(PWLS)(ρ(r))=

_(WLS)(ρ(r))+κ∥ρ(r)−ρ_(r)(r)|_(I) _(n) ²  (6)where I_(R) is the R-dimensional identity matrix. The regularizationparameter κ>0 ensures a balance between the closeness to the data andthe penalty term, which controls the deviation from a given referencevector ρ_(r)(r). The solution {circumflex over (ρ)}_(PWLS)(r) admits thefollowing closed-form expression:{circumflex over (ρ)}_(PWLS)(r)=ρ_(r)(r)+(S ^(H)(r)Ψ⁻¹ S(r)+κI _(R))⁻¹ S^(H)(r)Ψ⁻¹(d(r)−S(r)ρ_(r)(r))  (7).Note that the accuracy of the solution depends on the reference vectorand the choice of the regularization parameter κ.

Tikhonov regularization is known to introduce blurring in the image.

To overcome this limitation, “edge-preserving” penalty terms have beenproposed, which are applied in the image domain and make theregularization more efficient by limiting blurring effects andpreserving the image boundaries. However, the introduction of theseterms may lead to a non-differentiable optimization problem which is notalways easy to solve numerically.

As mentioned above, in [Chaâri et al. 2008] and [Chaâri et al. 2009] thepresent inventors have proposed a new regularization scheme; based onwavelet transforms (WT); which will be summarized here.

Indeed, in SENSE-based reconstructed images, well spatially-localizedartifacts appear as distorted curves with either very high or very lowintensity, and the WT has been recognized as a powerful tool thatenables a good space and frequency localization of useful information.

A general introduction to wavelets and wavelet transforms is provided by[Pesquet-Popescu, Pesquet].

In what follows, T stands for the WT operator. It corresponds to adiscrete decomposition onto a separable 2D M-band wavelet basisperformed over j_(max) resolution levels. The objective image ρ of sizeY×X can be viewed as an element of the Euclidean space

^(K) with K=Y×X endowed with the standard inner product

·|·

and norm ∥·∥. Let (e_(k))_(1≤k≤K) be the considered discrete waveletbasis of the space

^(K). The wavelet decomposition operator T is defined as the linearoperator:T:

^(K)→

^(K)ρ

(

ρ|e _(k)

)_(1≤k≤K)

The adjoint operator T* serving for reconstruction purposes is thendefined as the bijective linear operator:

${T^{*}\text{:}{\mathbb{C}}^{K}}->\left. {{\mathbb{C}}^{K}\left( \zeta_{k} \right)}_{1 \leq k \leq K}\mapsto{\sum\limits_{k = 1}^{K}{\zeta_{k}{e_{k}.}}} \right.$The resulting wavelet coefficient field of a target image function ρ isdefined by ζ=((ζ_(a,k))_(1≤k≤Kjmax), (ζ_(o,j,k))_(1≤j≤jmax, 1≤k≤K) _(j)) where K_(j)=KM^(−2j) is the number of wavelet coefficients in a givensubband at resolution j (by assuming that Y and X are multiple ofM^(jmax)) and the coefficients have been reindexed in such a way thatζ_(a,k) denotes an approximation coefficient at resolution level j_(max)and ζ_(o,j,k) denotes a detail coefficient at resolution level j andorientation o∈0={0, . . . , M−1}²\{(0,0)}.

In the dyadic case (M=2), there are three orientations corresponding tothe horizontal, vertical or diagonal directions. When an orthonormalwavelet basis is considered, the adjoint operator T* reduces to theinverse WT operator T⁻¹ and the operator norm ∥T∥ of T is equal to 1.

An estimate of the target image ρ is generated through thereconstruction wavelet operator T*. Let ζ be the unknown waveletcoefficients such that ρ=T*ζ. The aim is to build an estimate{circumflex over (ζ)} of the vector of coefficients ζ from theobservations d. This is based on a Bayesian approach relying on suitablepriors on the wavelet coefficients.

Given the observation model in Eq. (2.) and the assumptions regardingthe noise (i.i.d. circular Gaussian with zero-mean and between-coilcorrelation matrix Ψ), the likelihood function factorizes over pixelslying in the Y×X FOV:

$\begin{matrix}{{{p\left( d \middle| {T^{*}\zeta} \right)} = {{\prod\limits_{r \in {{\{{1,\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\ldots\mspace{14mu},X}\}}}}{p\left( {d(r)} \middle| {\rho(r)} \right)}} \propto {\prod\limits_{r \in {{\{{1,\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\ldots\mspace{14mu},X}\}}}}{\exp\left( {\mathcal{J}_{WLS}\left( {\rho(r)} \right)} \right)}}}}\mspace{20mu}{{{where}\mspace{14mu}\rho} = {T^{*}\zeta\mspace{14mu}{and}}}} & (8) \\{\mspace{20mu}{{\mathcal{J}_{L}(\rho)} = {\sum\limits_{r \in {{\{{1,\ldots\mspace{14mu},Y}\}} \times {\{{1,\ldots\mspace{14mu},X}\}}}}{\mathcal{J}_{WLS}\left( {\rho(r)} \right)}}}} & (9)\end{matrix}$

Let ƒ be the prior probability density function (pdf) of the image inthe wavelet domain. It will be assumed here that the real and imaginaryparts of the wavelet coefficients are independent. It will also beassumed that the real (resp. imaginary) parts of the waveletcoefficients are independent and identically distributed (i.i.d.) ineach subband. Their statistical characteristics may however vary betweentwo distinct subbands. Furthermore, by looking at the empiricaldistributions of the real and imaginary parts of the considered waveletcoefficients, the present inventors have noticed that their empiricalhistograms are well-fitted by a “Generalized Gauss-Laplace” (GGL)distribution, which presents a single mode and whose shape variesbetween the Gaussian and Laplacian densities. The corresponding pdfreads:

$\begin{matrix}{{\forall{\xi \in {{\mathbb{R}}.{f\left( {\xi\text{:}{\alpha \cdot \beta}} \right)}}}} = {\sqrt{\frac{\beta}{2\pi}}\frac{e^{- {({{\alpha{\xi }} + {\frac{\beta}{2}\zeta^{2}} + \frac{\sigma^{2}}{2\beta}})}}}{{erfc}\left( \frac{\alpha}{\sqrt{2\beta}} \right)}}} & (10)\end{matrix}$

where α∈

₊ and β∈

*₊ are hyper-parameters to be estimated. FIG. 1 illustrates theempirical histograms of real and imaginary parts of the horizontaldetail subband at the first resolution level using the dyadic (M=2)wavelet decomposition with Daubechies filters of length 8. This figureshows also that the adopted GGL distribution better fits the empiricalhistogram than a Generalized Gaussian (GG) pdf.

At the coarsest resolution level j_(max), the distributions of both thereal and imaginary parts of the approximation coefficients can bemodeled by a Gaussian distribution since they belong to a low frequencysubband.

Due to its familiarity and simplicity, the MAP (Maximum A Posteriori)estimator will be used for the estimation purpose, although differentchoices are possible. Based on the prior and the likelihood givenhereabove, the MAP estimator is computed by maximizing the fullposterior distribution or minimizing its negative log-likelihood:

$\begin{matrix}{{\hat{\zeta}}^{MAP} = {\underset{\zeta \in C^{K}}{\arg\mspace{11mu}\max}\left( {{\ln\;{p\left( d \middle| {T^{*}\zeta} \right)}} + {\ln\;{f(\zeta)}}} \right)}} & (11)\end{matrix}$or equivalently by minimizing the following criterion:

$\begin{matrix}\begin{matrix}{{\hat{\zeta}}^{MAP} = {\underset{\zeta \in {\mathbb{C}}^{K}}{\arg\mspace{11mu}\min}{\mathcal{J}_{WT}(\zeta)}}} \\{= {\underset{\zeta \in {\mathbb{C}}^{K}}{\arg\mspace{11mu}\min}\left( {{\mathcal{J}_{L}\left( {T^{*}\zeta} \right)} + {\mathcal{J}_{P}(\zeta)}} \right)}}\end{matrix} & (12)\end{matrix}$with:

$\begin{matrix}{\mspace{20mu}{{{{\mathcal{J}_{P}(\zeta)} = {{\sum\limits_{k = 1}^{K_{j_{{ma}\; x}}}{\Phi_{a}\left( \zeta_{a,k} \right)}} + {\sum\limits_{o \in O}{\sum\limits_{j = 1}^{j_{{{ma}\; x}\;}}{\sum\limits_{k = 1}^{K_{j}}{\Phi_{o,j}\left( \zeta_{o,j,k} \right)}}}}}}\mspace{20mu}{{\Phi_{a}\left( \zeta_{a,k} \right)} = {\frac{\left( {{{Re}\left( \zeta_{a,k} \right)} - \mu^{Re}} \right)^{2}}{\left( {\sqrt{2}\sigma^{Re}} \right)^{2}} + \frac{\left( {{{Im}\left( \zeta_{a,k} \right)} - \mu^{Im}} \right)^{2}}{\left( {\sqrt{2}\sigma^{Im}} \right)^{2}}}}},{{\Phi_{o,j}\left( \zeta_{o,j,k} \right)} = {{{\alpha_{o,j}^{Re}{{{Re}\left( \zeta_{o,j,k} \right)}}} + {\frac{\beta_{o,j}^{Re}}{2}{{{Re}\left( \zeta_{o,j,k} \right)}}^{2}} + {\alpha_{o,j}^{Im}{{{Im}\left( \zeta_{o,j,k} \right)}}} + {\frac{\beta_{o,j}^{Im}}{2}{{{Im}\left( \zeta_{{o,j,k}\;} \right)}}^{2}}}..}}}} & (13)\end{matrix}$

Hereabove, Re(·) and Im(·) (or ·^(Re) and ·^(Im)) stand for the real andimaginary parts, respectively. The prior parameters (hyperparameters)α_(o,j)=(α_(o,j) ^(Re), α_(o,j) ^(Im))∈(

*₊)², β_(o,j)=(β_(o,j) ^(Re), β_(o,j) ^(Im))∈(

*₊)², μ_(o,j)(μ_(o,j) ^(Re), μ_(o,j) ^(Im))∈(

)² and σ_(o,j)=(σ_(o,j) ^(Re), σ_(o,j) ^(Im))∈(

₊)² are unknown and need to be estimated. The estimation ofhyperparameters is an important part of the invention and it will bediscussed extensively later.

While J_(L) is differentiable with a Lipschitz-continuous gradient,J_(P) is not differentiable. Therefore, although the penalty functionJ_(WT) is convex, the optimization procedure cannot rely on conventionalconvex optimization techniques like the pseudo-conjugate gradient.

Optimization can be performed using a generalized form of the iterativeoptimization procedure developed in [Daubechies et al., 2004; Chaux etal., 2007], which is based on the Forward-Backward (FB) algorithm and onthe concept of proximity operator, generalized to the case of functionsof a complex variable.

For the functionϕ:

^(K)→]−∞,+∞]x

ϕ ^(Re)(Re(x))+ϕ^(Im)(x))  (14)

where ϕ^(Re) and ϕ^(Im) are functions in Γ₀(R^(K)) and Re(x) (resp.Im(x)) is the vector of the real parts (resp. imaginary parts) of thecomponents of x∈

^(K) the proximity operator is defined asprox_(ϕ):

^(K)→

^(K) x

prox_(ϕRe)(Re(x))+ιprox_(ϕIm)(Im(x)).

By extending the algorithm in [Chaux et al., 2007] to the complex case,a minimizer of J_(WT) can be iteratively computed according to thealgorithm below (Algorithm 1), where the gradient of J_(L) is firstcalculated, and then the frame coefficients are updated. It can benoticed that λ_(n) and γ_(n) correspond to relaxation and step-sizeparameters, respectively.

Algorithm 1: Let (γ_(n))_(n>0) and (λ_(n))_(n>0) sequences of positivereals.   1: Set n = 0 and ε ≥ 0. Initialize ζ^((n)) and set

 ^((n)) =

 _(WT)(ζ^((n))).   2: repeat   3:  Reconstruct the image by settingρ^((n)) = T * ζ^((n)).   4:  Compute the image u^((n)) such that:  ∀r ∈{1 , . . . , Y/R} × {1 , . . . , X},  u^((n))(r) = 2S^(H)(r)Ψ⁻¹(S(r)ρ^((n))(r) − d(r)),  where the vector u^((n)) (r) is defined fromu^((n)) in the same way as ρ(r) is defined  from  ρ (see Eq. (2.17)).  5:  Determine the wavelet coefficients ν^((n)) = Tu^((n)) = (ν_(α),(ν_(o,j)

_(1≤j≤j) _(max) ) of u^((n)).   6:  Update the approximationcoefficients of the reconstructed image ρ^((n+1)):  ∀k ∈ {1, . . . ,K_(j) _(max) }, ζ_(α,k) ^((n+1)) = ζ_(α,k) ^((n)) +λ_(n)(prox_(γnϕα)(ζ_(α,k) ^((n)) − γ_(n)ν_(α,k) ^((n))) − ζ_(α,k)^((n))).   7:  Update the detail coefficients of the reconstructed imageρ^((n+1)):   ∀o ∈  

, ∀j ∈ {1, . . . , j_(max)}, ∀k ∈ {1, . . . , K_(j)},    ζ_(o,j,k)^((n+1)) = ζ_(o,j,k) ^((n)) + λ_(n)(prox_(γnϕo,j)(ζ_(o,j,k) ^((n)) −γ_(n)ν_(o,j,k) ^((n))) − ζ_(o,j,k) ^((n))).   8:  Compute

 ^((n+1)) =

 _(WT)(ζ^((n+1))).   9:  n ← n + 1  10: until | 

 ^((n)) −

 ^((n−1))| ≤ ε 

 ^((n−1))  11: return ρ^((n)) = T * ζ^((n))

For every r⊂{1, . . . , Y}×{1, . . . , X}, let θ_(r)≥0 be the maximumeigenvalue of the Hermitian positive semi-definite matrix S^(H)(r)Ψ⁻¹S(r) and let θ=max_(r∈{1, . . . , Y}×{1, . . . , X})θ_(r)>0. Toguarantee the convergence of Algorithm 1, the step-size and relaxationparameters are subject to the following conditions:

${{(i)\mspace{14mu}\inf_{n > 0}\gamma_{n}} > {0\mspace{14mu}{and}\mspace{14mu}\sup_{n > 0}\gamma_{n}} < \frac{1}{\theta{T}^{2}}},{{({ii})\mspace{14mu}\inf_{n > 0}\lambda_{n}} > {0\mspace{14mu}{and}\mspace{14mu}\sup_{n > 0}\lambda_{n}} \leq 1.}$

As discussed in [Chaâri et al. 2009], the optimization problem (11),(12) can be modified to include additional (preferably convex)constraints, e.g. upper and lower bounds of the image intensity value.

The present inventors have found that using wavelet regularization, asdiscussed above, allows preserving the image details while smoothingreconstruction artifacts but may introduce some irregularity inhomogeneous area of the image. On the other hand, other (nonwavelet-based) regularization schemes are known to be adapted toregularize smooth regions, but at the expense of an oversmoothing ofimage details. This is particular the case of “Total Variation” TVregularization, described inter alia by [Raj et al., 2007] and [Block etal., 2007].

A first improvement of the method of [Chaâri et al. 2008] and [Chaâri etal. 2009], constituting an aspect of the present invention, consists incombining WT and TV in a joint regularization framework in order to takeadvantage of their different properties, allowing them to alleviate thedrawbacks of each other. A joint Wavelet Transform-Total Variation(WTTV) regularization can reduce to the optimization of the followingpenalized criterion:

$\begin{matrix}\begin{matrix}{{\hat{\rho}}_{{WT} - {TV}} = {T^{*}\underset{\zeta}{\arg{\;\;}\min}{\mathcal{J}_{{WT} - {TV}}(\zeta)}}} \\{= {{T^{*}\underset{\zeta\;}{\arg\mspace{11mu}\min}{\mathcal{J}_{L}\left( {T^{*}\zeta} \right)}} + {\kappa_{1}{\mathcal{J}_{P}(\zeta)}} + {\kappa_{2}{{T^{*}\zeta}}_{TV}}}}\end{matrix} & (15)\end{matrix}$where J_(L) and J_(P) have bee, defined above, κ₁>0 and κ₂>0 areregularization parameters, and T* is the wavelet adjoint operator. Thediscrete total variation norm of an image

P _(E) cYxX norm of an image ρ∈

^(Y×X) is given by:

$\begin{matrix}{{\rho }_{TV} = {\sum\limits_{y = 1}^{Y}\;{\sum\limits_{x = 1}^{X}\;\sqrt{{{\left( {\nabla_{1}\rho} \right)\left( {y,x} \right)}}^{2} + {{\left( {\nabla_{1}\rho^{T}} \right)^{T}\left( {y,x} \right)}}^{2}}}}} & (16)\end{matrix}$

where for every ρ∈

^(Y×X), ∇₁ is the horizontally smoothed gradient operator defined by

$\begin{matrix}{{\nabla_{1}(\rho)} = {{\frac{1}{2}\left( {{\rho\left( {{y + 1},{x + 1}} \right)} - {\rho\left( {y,{x + 1}} \right)} + {\rho\left( {{y + 1},x} \right)} - {\rho\left( {y,x} \right)}} \right)}:_{{\leq y \leq Y},{1 \leq x \leq X}}}} & (17)\end{matrix}$

For the sake of presentation, it has been assumed that the image ρ isperiodic or equivalently that the image boundaries are toroidal.

Minimizing the optimality criterion (15) is much more difficult thanminimizing the space-only criteria given by (11)-(13), because more thantwo terms are involved, while the Forward-Backward (FB) algorithmdiscussed above only applies to the minimization of criteria comprisingtwo terms. Iterative minimization of a non-differentiable convexfunction comprising more than two terms can be performed using theso-called Parallel ProXimal Algorithm (PPXA) [Combettes and Pesquet,2008], which also requires calculating the proximity operator of each ofthree involved terms.

The difficulty stems from the calculation of the proximity operator of∥·∥_(TV)oT*. To circumvent this difficulty, the TV penalization inequation (15) can be split into four terms as [Combettes and Pesquet,2008]. The TV penalty term therefore reads:

${\forall{\rho \in {\mathbb{C}}^{Y \times X}}},{{\rho }_{TV} = {\sum\limits_{i = 0}^{3}\;{{tv}_{i}(\rho)}}}$where for every (q, r)∈{0,1}²,

${{tv}_{{2q} + r}(\rho)} = {\sum\limits_{y = 1}^{Y/2}\;{\sum\limits_{x = 1}^{X/2}\;\sqrt{\begin{matrix}{{{\left( {\nabla_{1}\rho} \right)\left( {{{2y} + q},{{2x} + r}} \right)}}^{2} +} \\{{\left( {\nabla_{1}\rho^{T}} \right)^{T}\left( {{{2y} + q},{{2x} + r}} \right)}}^{2}\end{matrix}}}}$For every q and r in {0, 1}, let ↓^(q,r) be the decimation operatordefined by↓_(q,r):

^(2Y×2X)→

^(Y×X)ν=ν_(y,x))_(1≤y≤2Y,1≤x≤2,X)

(ν_(2y+q,2r+r))_(1≤y≤2Y,1≤x≤2X).

and U_(q+2r) be the following operator:

U_(q + 2r) : ℂ^(Y × X) → ℂ^(Y × X)$\left. \rho\mapsto\left. \downarrow{}_{q,r}\begin{bmatrix}{\nabla_{0}\rho} & {\nabla_{1}\rho} \\\left( {\nabla_{1}\rho^{T}} \right)^{T} & {\nabla_{2}\rho}\end{bmatrix} \right. \right.$where for every ρ∈

^(Y×X), the operators ∇₀ and ∇₂ are defined by

${\nabla_{0}(\rho)} = {\frac{1}{2}\left( {{\rho\left( {{y + 1},{x + 1}} \right)} + {\rho\left( {y,{x + 1}} \right)} + {\rho\left( {{y + 1},x} \right)} + {\rho\left( {y,x} \right)}} \right)_{{1 \leq y \leq Y},{1 \leq x \leq X}}}$  and${\nabla_{2}(\rho)} = {\frac{1}{2}\left( {{\rho\left( {{y + 1},{x + 1}} \right)} + {\rho\left( {y,{x + 1}} \right)} - {\rho\left( {{y + 1},x} \right)} + {\rho\left( {y,x} \right)}} \right)_{{1 \leq y \leq Y},{1 \leq x \leq X}}}$

Let also h be the function defined on

^(Y×X) by:

${h(\rho)} = {\sum\limits_{y = 1}^{Y/2}\;{\sum\limits_{x = 1}^{X/2}\;\sqrt{{{\left( {\nabla_{1}\rho} \right)\left( {y,{x + {X/2}}} \right)}}^{2} + {{\left( {\nabla_{1}\rho^{T}} \right)^{T}\left( {{y + {Y/2}},x} \right)}}^{2}}}}$

It turns out from the equations above that ∀i∈{0, 1, 2, 3},tν_(i)=hoU_(i). Consequently, the optimization problem in equation (15)can be rewritten as:

$\begin{matrix}\begin{matrix}{{\hat{\rho}}_{{WT} - {TV}} = {T^{*}{\underset{\zeta}{\arg\mspace{11mu}\min}\left\lbrack {{\mathcal{J}_{L}\left( {T^{*}\zeta} \right)} + {\kappa_{1}{\mathcal{J}_{P}(\zeta)}} + {\kappa_{2}{\sum\limits_{i = 0}^{3}\;{{tv}_{i}\left( {T^{*}\zeta} \right)}}}} \right\rbrack}}} \\{= {T^{*}{\underset{\zeta}{\arg\mspace{11mu}\min}\left\lbrack {{\mathcal{J}_{L}\left( {T^{*}\zeta} \right)} + {\kappa_{1}{\mathcal{J}_{P}(\zeta)}} + {\kappa_{2}{\sum\limits_{i = 0}^{3}\;{h\left( {U_{i}\left( {T^{*}\zeta} \right)} \right)}}}} \right\rbrack}}}\end{matrix} & (18)\end{matrix}$

The regularized WT-TV reconstruction approach is summarized in Algorithm2 where the PPXA algorithm [Combettes and Pesquet, 2008] is used tominimize the optimality criterion in equation (18) which is made up ofJ=6 convex functions.

Algorithm 2 Set γ ϵ]0, + ∞[, n = 0, (ω_(i))_(1≤i≤6) ∈ [0,1]⁶ such thatΣ_(i=0) ⁶ ωi = 1, (ζ_(i) ^((n)))1 ≤ i ≤ 6 ∈ (

^(K))⁶ where ζ_(i) ^((n)) = ((ζ_(α) ^(i,(n))), (ζ_(o,j) ^(i,(n))))o ϵ

_(,)1 ≤ j ≤ j _(max) ) for every i ∈ {1, . . . , 6}. Set also ε ≥ 0_(χ)initialize ζ^((n)) = Σ_(i=1) ⁶ ω_(i) ζ_(i) ^((n)) and

^((n)) = 0.  1: repeat  2: Compute the image u^((n)) such that: ∀r ∈ {1,. . . , Y / R} × {1, . . . , X},${u^{(n)}(r)} = {{\left( {I_{R} + {\frac{2\gamma}{\omega_{1}}{S^{H}(r)}\Psi^{- 1}{S(r)}}} \right)^{- 1}\left( {{\rho_{1}^{(n)}(r)} + {\frac{2\gamma}{\omega_{1}}{S^{H}(r)}\Psi^{- 1}{d(r)}}} \right)\mspace{14mu}{where}\mspace{14mu}\rho_{1}^{(n)}} = {T^{*}{\zeta_{1}^{(n)}.}}}$ 3: Compute the wavelet coefficients p₁ ^((n)) = Tu^((n)) of u^((n)). 4: Calculate p₂ ^((n)) = {prox_(γK1Φ) _(α/ω2) (ζ_(α) ^(2,(n))),(prox_(γΦ) _(o,j/ω2) (ζ_(o,j) ^(2,(n))))_(o∈)

_(,1≤j≤j) _(max) ).  5: For every i ∈ {0,1,2,3}, calculate p_(i+3)^((n)) = prox_(γK2tv) _(i/ωi) (ζ_(i+3) ^((n))).  6: Set P^((n)) =Σ_(i=1) ⁶ ω_(i)P_(i) ^((n)).  7: Set λ_(n) ∈ [0,2].  8: for i = 1 to 6do  9: ζ_(i) ^((n)) = ζ_(i) ^((n)) + λ_(n)(2P^((n)) − ζ^((n)) − P_(i)^((n))). 10: end for 11: ζ^((n+1)) = ζ^((n)) + λ_(n)(P^((n)) − ζ^((n))).12: Compute

^((n+1)) =

_(WT−TV)(ζ^((n+1))). 13: n ← n + 1. 14: until |

^((n)) −

^((n−1))| ≤ ε

^((n−1)) 15: return ρ^((n)) = T * ζ^((n))

In the method of [Chaâri et al. 2008] and [Chaâri et al. 2009],three-dimensional (3D) imaging is performed by stacking regularizedtwo-dimensional (2D) images. A further improvement of the method,constituting another aspect of the present invention, consists inperforming direct, regularized 3D image reconstruction. In the 3D case,the image reconstruction problem can still be written:d(r)=S(r)ρ(r)+n(r)  (19)like in equation (12); however, r=(y, x, z) is now the three-dimensionalspatial position, z∈{1, . . . , Z} being the position along the thirddirection (slice selection). Moreover, the penalty term J_(P) willdepend on the distribution of the wavelet coefficient of a 3D dyadicwavelet transform.

A further improvement of the method, constituting another aspect of thepresent invention, consists in performing four-dimension (4D),spatio-temporal regularization in dynamic MRI. This embodiment of theinvention can be applied e.g. to fMRI of the brain, wherein the wholebrain volume has to be acquired several times, yielding a 4D dataset. Inconventional fMRI, the 3D images are supposed independent although theybelong to the same fMRI session. However, in practice, the 3D temporalimages are somehow dependent since they belong to the same fMRI sessioninvolving the same experimental paradigm. The BOLD (Blood Oxygen LevelDependent) time-series and the acquisition noise are in fact correlatedin time in a brain fMRI session. For this reason, taking into accounttemporal dependencies between 3D images helps to increase the SNRthrough the acquired images, and therefore enhances the reliability ofthe statistical analysis in fMRI. However, since in dynamic MRI theimaged object geometry generally changes during the acquisition, joiningthe reconstruction process to the temporal regularization is verydifficult.

To deal with a 4D reconstruction of a time series of N_(r)three-dimensional images (N_(r) is usually even), the observation modelof equations (12) and (19) will be rewritten as follows:d _(t)(r)=S(r)ρ_(t)(r)+n _(t)(r)  (20)where t∈{1, . . . , N_(r)} is the acquisition time. Using a dyadic 3Dwavelet operator T, the coefficients will be reindexed so thatζ^(t)=(ζ_(α) ^(t), (ζ_(o,j) ^(t))_(o∈)

_(,1≤j≤j) _(max) ) with o∈

={0,1}³\{(0,0,0)}.

Accounting for an additional temporal l_(p) regularization term,reconstruction of the 4D “volume” (i.e. of the time series of 3D imagesor volumes) is performed through the minimization of the followingoptimality criterion:

$\begin{matrix}\begin{matrix}{\hat{\zeta} = {\arg{\;\;}{\min\limits_{\zeta}{\mathcal{J}_{ST}(\zeta)}}}} \\{= {\arg\mspace{11mu}{\min\limits_{\zeta}{\sum\limits_{t = 1}^{N_{r}}\;\sum\limits_{r \in {{\{{1,\mspace{11mu}\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},X}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},Z}\}}}}}}}} \\{{{{d_{t}(r)} - {{S(r)}\left( {T^{*}\zeta^{t}} \right)(r)}}}_{\Psi^{- 1}}^{2} +} \\{{\sum\limits_{t = 1}^{N_{r}}\;{\mathcal{J}_{P}\left( \zeta^{t} \right)}} + {\kappa{\sum\limits_{t = 2}^{N_{r}}\;{{{T^{*}\zeta^{t}} - {T^{*}\zeta^{t - 1}}}}_{p}^{p}}}}\end{matrix} & (21)\end{matrix}$where ζ=(ζ¹, ζ², . . . , ζ^(Nr))T, κ>0 is a regularization parameter andJ_(P) is defined as in (13).

In equation (21),

$\sum\limits_{t = 1}^{N_{r}}\;{\sum\limits_{r \in {{\{{1,\mspace{11mu}\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},X}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},Z}\}}}}\;{{{d_{t}(r)} - {{S(r)}\left( {T^{*}\zeta^{t}} \right)(r)}}}_{\Psi^{- 1}}^{2}}$is the “error term” representative of a likelihood of each reconstructedimage, given corresponding acquired elementary images;

$\sum\limits_{t - 2}^{N_{r}}\;{{{{T^{*}\zeta^{t}} - {T^{*}\zeta^{t - 1}}}}_{p}^{p}.}$is a temporal penalty term, representative of a pixel-by-pixel orvoxel-by-voxel difference between consecutive images of the series and

$\sum\limits_{t = 1}^{N_{r}}\;{\mathcal{J}_{P}\left( \zeta^{t} \right)}$is a wavelet penalty term used for spatial regularization. It will beunderstood that a different (e.g. non wavelet-based) spatialregularization term could also be used.

The adjoint wavelet operator T* is then applied to each component ζ^(t)of ζ to obtain the reconstructed 3D image ρ^(t) at the acquisition timet by taking into account the time dependencies with the other acquiredimages.

It will be noted that the temporal penalty term is not based onwavelets, but is representative of a pixel-by-pixel or voxel-by-voxeldifference between consecutive images of the series. More precisely, inequation (21) the temporal penalty term is based on an L_(p) norm; thepresent inventors have found that p has to be greater or equal to one,and that it should preferably satisfy: 1≤p<1.5. Other form of thetemporal penalty term could be used, in particular based on (preferablyconvex) edge-preserving functions of a pixel-by-pixel or voxel-by-voxeldifference between consecutive images of the series.

Minimizing the optimality criterion (21) is much more difficult thanminimizing the space-only criteria given by (11)-(13), because more thantwo terms are involved, while the Forward-Backward (FB) algorithmdiscussed above only applies to the minimization of criteria comprisingtwo terms. Iterative minimization of a non-differentiable convexfunction comprising more than two terms can be performed using thealready-cited Parallel ProXimal Algorithm (PPXA) [Combettes and Pesquet,2008], which also requires calculating the proximity operator of each ofthree involved terms. This task is quite simple for the two first termsof equation (21) since they are separable with respect to the timevariable t and the spatial position r. However, this is not the case forthe time penalization term (third term in Eq. (21), which makesnon-trivial the calculation of the corresponding proximity operator. Itis then proposed to rewrite the optimality criteria in J_(ST) bydecomposing the time penalization into two terms which are separablewith respect to the time variable t (a given acquisition time t iseither involved in the first or in the second term), and for which theproximity operators are easy to calculate. More specifically, thetemporal penalty term J_(T) is expressed as the sum of a first partialtemporal penalty term J_(T) ¹ and a second partial temporal penalty termJ_(T) ², wherein: the first partial temporal penalty term isrepresentative of pixel-by-pixel or voxel-by-voxel differences betweeneach even-numbered image of the series and a preceding odd-numberedimage; and the second partial temporal penalty term is representative ofpixel-by-pixel or voxel-by-voxel differences between each odd-numberedimage of the series and a preceding even-numbered image:

$\begin{matrix}{{{\mathcal{J}_{ST}(\zeta)} = {{\sum\limits_{t = 1}^{N_{r}}\;{\sum\limits_{r \in {{\{{1,\mspace{11mu}\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},X}\}} \times {\{{1,\mspace{11mu}\ldots\mspace{14mu},Z}\}}}}{{{d_{t}(r)} - {{S(r)}\left( {T^{*}\zeta^{\ell}} \right)(r)}}}_{\Psi^{- 1}}^{2}}} + {\sum\limits_{t = 1}^{N_{r}}\;{\mathcal{J}_{P}\left( \zeta^{t} \right)}} + {\mathcal{J}_{T}^{1}(\zeta)} + {\mathcal{J}_{T}^{2}(\zeta)}}}\mspace{20mu}{{\mathcal{J}_{T}^{1}(\zeta)} = {\kappa{\sum\limits_{t = 1}^{N_{r}/2}\;{{{T^{*}\zeta^{2t}} - {T^{*}\zeta^{{2t} - 1}}}}_{p}^{p}}}}\mspace{20mu}{{\mathcal{J}_{T}^{2}(\zeta)} = {\kappa{\sum\limits_{t = 1}^{{N_{r}/2} - 1}\;{{{T^{*}\zeta^{{2t} + 1}} - {T^{*}\zeta^{2t}}}}_{p}^{p}}}}} & (22)\end{matrix}$

Since J_(T) ¹ and J_(T) ² are separable with respect to the timevariable t, the corresponding proximity operator can easily becalculated based on the proximity operator of each of the involvedterms. Let us consider the following functionΨ:

^(K)×

^(K)→

(ζ^(t),ζ^(t−1))

κ∥T*ζ ^(t) −T*ζ ^(t−1)∥_(p) ^(p) =ϕoH((ζ^(t),ζ^(t−1)))where ϕ(·)=κ∥T*·∥_(p) ^(p) and H is a linear operator defined byH:

^(K)×

^(K)→

^(K)(a,b)

a−b.

Its associated adjoint operator H* is therefore given byH*:

^(K)→

^(K)×

^(K) a

(a−a).

The proximity operator of Φ is then given by:

$\begin{matrix}{{prox}_{\Phi} = {{prox}_{\phi\; \circ H} = {{Id} + {\frac{1}{2}{{H^{*}\left( {{prox}_{2\phi} - {Id}} \right)} \circ H}}}}} & (23)\end{matrix}$

The resulting algorithm (Algorithm 3) for the minimization of thespace-time optimality criterion J_(ST) is given below.

Algorithm 3 Set γ ϵ]0, + ∞[, n = 0, (ω_(i))_(1≤i≤4) ∈ [0,1]⁴ such thatΣ_(i=0) ⁴ ω_(i) = 1, (ζ_(i) ^((n)))_(1≤i≤4) ∈ (

^(K×N) ^(r) )⁴ where ζ_(i) ^((n)) = (ζ_(i) ^(1,(n)), ζ_(i) ^(2,(n)), . .. , ζ_(i) ^(Nr,(n))) and ζ_(i) ^(t,(n)) = ((ζ_(i,a) ^(t,(n))),((ζ_(i,o,j) ^(t,(n))))o ϵ

_(,)1 ≤ j ≤ j_(max)) for every i ∈ {1, . . . , 4} and t ∈ {1, . . . ,N_(r)} Set also ε ≥ 0 and initialize ζ^((n)) = Σ_(i=1) ⁴ ω_(i)ζ_(i)^((n)) and

^((n)) = 0. 1: repeat 2: Set p₄ ^(1,(n)) = ζ₄ ^(1,(n)). 3: for t = 1 toN_(r) do 4: Calculate the image u_(t) ^((n)) such that: ∀r ∈ {1, . . . ,Y/R} × {1, . . . , X} × {1, . . . , Z},${u_{t}^{(n)}(r)} = {{\left( {I_{R} + {\frac{2\gamma}{\omega_{1}}{S^{H}(r)}\Psi^{- 1}{S(r)}}} \right)^{- 1}\left( {{\rho_{t,1}^{(n)}(r)} + {\frac{2\gamma}{\omega_{1}}{S^{H}(r)}\Psi^{- 1}{d_{t}(r)}}} \right)\mspace{14mu}{where}\mspace{14mu}\rho_{t,1}^{(n)}} = {T^{*}{\zeta_{1}^{t,{(n)}}.}}}$5: Compute the wavelet coefficients p₁ ^(t,(n)) = Tu_(t) ^((n)). 6:Compute p₂ ^(t,(n)) = (prox_(γΦ) _(α/ω2) (ζ_(2,α) ^(t,(n))),(prox_(γΦ)_(o,j/ω2) (ζ_(2,o,j) ^(t,(n)))) o ∈

_(,)1 ≤ j ≤ j_(max). 7: if t is even then 8: calculate (p₃ ^(t,(n)), p₃^(t−1,(n))) = prox_(γKΦ/ω3)((ζ^(t,(n)), ζ^(t−1(n)))) 9: else if t is oddand t > 1 then 10: calculate (p₄ ^(t,(n)), p₄ ^(t−1,(n))) =prox_(γKΦ/ω4)((ζ^(t,(n)), ζ^(t−1,(n)))). 11: end if 12: if t > 1 then13: Set P^(t−1,(n)) = Σ_(i=1) ⁴ ω_(i)p_(i) ^(t−1,(n)). 14: end if 15:end for 16: Set p₄ ^(Nr,(n)) = ζ₄ ^(Nr,(n)). 17: Set P^(Nr,(n)) =Σ_(i=1) ⁴ ω_(i)p_(i) ^(Nr,(n)). 18: Set p₁ ^((n)) = (p₁ ^(1,(n)), p₁^(2,(n)), . . . , p₁ ^(Nr,(n))), p₂ ^((n)) = (p₂ ^(1,(n)), p₂ ^(2,(n)),. . . , p₂ ^(Nr,(n))), p₃ ^((n)) = (p₃ ^(1,(n)), p₃ ^(2,(n)), . . . , p₃^(Nr,(n))), p₄ ^((n)) = (p₄ ^(1,(n)), p₄ ^(2,(n)), . . . , p₄ ^(Nr,(n)))and P^((n)) = (P^(1,(n)), P^(2,(n)), . . . P^(Nr,(n))). 19: Set λ_(n) ∈[0,2]. 20: for i = 1 to 4 do 21: ζ_(i) ^((n)) = ζ_(i) ^((n)) +λ_(n)(2P^((n)) − ζ^((n)) − p_(i) ^((n))). 22: end for 23: ζ^((n+1)) =ζ^((n)) + λ_(n)(P^((n)) − ζ^((n))). 24: Compute

^((n+1)) =

_(ST)(ζ^((n+1))). 25: n ← n + 1. 26: until |

^((n)) −

^((n−1))| ≤ ε

^((n−1)). 27: Set {circumflex over (ζ)} = ζ^((n)).

It will be understood that space-time regularization can also be appliedto the case of a time series of bi-dimensional images, even if this isless common. Moreover, the spatial penalty term may also include a totalvariation component and/or constraints, as in the case of purely spatial2D or 3D regularization.

It is worth noticing that in a 4D fMRI dataset, and since the imagedbody volume is acquired many times, said body (i.e. a brain) mayslightly move between successive scans. The inherent motion artifactsmay be dramatic for the fMRI analysis. In this context, and since astandard fMRI study involves a motion correction step through applying arigid transformation, another extension of the estimation method isproposed. This method accounts for the eventual motion artifacts byinvolving the rigid transformation parameters Sr in the estimation step.The transformation is first applied to the standard SENSE data whichwill be used for the estimation step. An extension of the abovedescribed methods is also proposed performing the estimation step not onthe standard SENSE reconstructed data, but on a quadraticallyregularized (QR) version of SENSE using Tikhonov regularization with asmoothing matrix (i.e. QR-SENSE). The QR-SENSE appears more efficientfor improving the robustness of the estimation step than introducing anadditional filtering step on the SENSE solution and remainscomputationally efficient since it admits a closed form expression. Inorder to further consider the motion artifacts occurring in fMRI, theproposed 4D reconstruction method (Algorithm 3) can also be extended toaccount for the rigid transformation parameters during thereconstruction process. The optimality criterion in equation (22) istherefore rewritten as follows:

$\begin{matrix}{{{\mathcal{J}_{ST}(\zeta)} = {{\sum\limits_{t = 1}^{N_{r}}\;{\sum\limits_{r \in {{\{{1,\;\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\;\ldots\mspace{14mu},X}\}} \times {\{{1,\;\ldots\mspace{14mu},Z}\}}}}\;{{{d_{t}\left( {r + {\delta\; r}} \right)} - {{S\left( {r - {\delta\; r}} \right)}\left( {T^{*}\zeta^{t}} \right)\left( {r + {\delta\; r}} \right)}}}_{\Psi^{- 1}}^{2}}} + {\sum\limits_{t = 1}^{N_{r}}\;{\mathcal{J}_{P}\left( \zeta_{c}^{t} \right)}} + {\mathcal{J}_{T}^{1}\left( \zeta_{c} \right)} + {\mathcal{J}_{T}^{2}\left( \zeta_{c} \right)}}},} & \left( {22{bis}} \right)\end{matrix}$

where ζ_(c) ^(t) is derived from ζ^(t) after applying the motioncorrection transformation and assuming that

${\mathcal{J}_{T}^{1}\left( \zeta_{c} \right)} = {\sum\limits_{t = 1}^{N_{r}/2}\;{\sum\limits_{r \in {{\{{1,\;\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\;\ldots\mspace{14mu},X}\}} \times {\{{1,\;\ldots\mspace{14mu},Z}\}}}}\;{\kappa_{r}{{{\left( {T^{*}\zeta_{c}^{2t}} \right)(r)} - {\left( {T^{*}\zeta_{c}^{{2t} - 1}} \right)(r)}}}^{p}}}}$  and${\mathcal{J}_{T}^{2}(\zeta)} = {\sum\limits_{t = 1}^{{N_{r}/2} - 1}\;{\sum\limits_{r \in {{\{{1,\;\ldots\mspace{14mu},{Y/R}}\}} \times {\{{1,\;\ldots\mspace{14mu},X}\}} \times {\{{1,\;\ldots\mspace{14mu},Z}\}}}}\;{\kappa_{r}{{{\left( {T^{*}\zeta_{c}^{{2t} + 1}} \right)(r)} - {\left( {T^{*}\zeta_{c}^{2t}} \right)(r)}}}^{p}}}}$Since

_(T) ¹ (resp.

_(T) ²) is separable with respect to the time variable t, its proximityoperator can easily be calculated based on the proximity operator ofeach of the involved terms in the sums above.

The rigid transformation parameters δr can be found by minimizing

_(ST)(ζ), expressed by equation (22 bis), also with respect to theseparameters. Alternatively, said rigid transformation parameters can beestimated separately, on the basis of a conventional SENSEreconstruction.

All the regularization methods which have been described here(wavelet-based 2D spatial regularization, wavelet-based 3D spatialregularization, hybrid wavelet-based/total variation 2D and 3D spatialregularization, 3D and 4D space-time regularization) involve numericalparameters which have to be set, either manually or automatically usinga suitable estimation algorithm. These parameters include the priorparameters, of the statistical distribution of the wavelet coefficients,or “hyperparameters”: α_(o,j)=(α_(o,j) ^(Re), α_(o,j) ^(Im))∈(

*₊)², β_(o,j)=(β_(o,j) ^(Re), β_(o,j) ^(Im))∈(

*₊)², μ_(o,j)=(μ_(o,j) ^(Re), μ_(o,j) ^(Im))∈(

)² and σ_(o,j)=(σ_(o,j) ^(Re), σ_(o,j) ^(Im))∈(

₊)². They can also include the temporal regularization parameter κ, theso-called “shape parameter”, i.e. the value of “p” defining the L_(p)norm ∥·∥_(p) used in temporal regularization (see equation 21), and theparameters κ₁, κ₂ determining the relative weight of the wavelet andtotal variation penalty terms in hybrid regularization.

The already-cited documents [Chaâri et al. 2008] and [Chaâri et al.2009] do not teach how to determine the relevant parameters. A furtherimprovement of the method described by these documents, constitutinganother aspect of the present invention, consists in providing automaticcalibration algorithms, allowing a fully or partly non-supervisedwavelet-based regularization. Two different frameworks are considered toachieve this goal depending on the type of decomposition. In case oforthonormal wavelet decomposition, the hyperparameters can be accuratelyestimated using a maximum likelihood procedure (see below). In contrast,when resorting to redundant frame decomposition, a stochastic samplingprocedure seems more efficient. In the latter case, the hyperparametersare estimated using the Minimum Mean Square Error (MMSE) or equivalentlythe posterior mean estimator (see below). The posterior mean estimatorcan also be used with an orthonormal wavelet decomposition, as this is aspecial case of frame.

First of all, the determination of the hyper-parameters defining theprior statistical distribution of the wavelet coefficients used for 2Dor 3D spatial regularization will be considered.

A first possibility for determining the hyper-parameters Θ=(μ, σ,(α_(o,j), β_(o,j))o∈

, 1≤j≤j_(max)) is to maximize the integrated likelihood according to thefollowing equation:

$\begin{matrix}{\hat{\Theta} = {{\underset{\Theta}{\arg\mspace{11mu}\max}{p\left( {d;\Theta} \right)}} = {\underset{\Theta}{\arg\mspace{11mu}\max}{\int{{p\left( d \middle| {T^{*}\zeta} \right)}{f\left( {\zeta;\Theta} \right)}d\;{\zeta.}}}}}} & (24)\end{matrix}$

Maximizing Eq. (24) is a missing data problem since ζ is unknown. Itrequires integrating out the sought image decomposition ζ and iteratingbetween image reconstruction and hyper-parameter estimation using theintensive EM algorithm described by [Dempster et al., 1977]. In order toalleviate the computational burden, it is advantageous to proceeddifferently by assuming that a reference or “auxiliary” full FOV image ρis available, and so is its wavelet decomposition ζ=(T*)⁻¹ ρ. Inpractice, the auxiliary image ρ can be obtained using 1D-SENSEreconstruction at the same R value, either unregularized orquadratically-regularized in a conventional way (e.g. using Tikhonovregularization) in the image space or in the k-space.

Then, the maximum likelihood (ML) estimation procedure consists ofassuming that this auxiliary image is a realization of the full priordistribution and thus in fitting Θ directly on it:

$\hat{\Theta} = {\underset{\Theta}{\arg\mspace{11mu}\max}\;{f\left( {\zeta;\Theta} \right)}}$In statistics, this solution is referred to as the complete data maximumlikelihood as opposed to the above mentioned missing data maximumlikelihood estimator. This procedure can be decomposed in twoindependent steps, the first one involving the setting of the Gaussianprior parameters (μ,σ) attached to the approximation coefficients ζ _(α)and the second one being related to the estimation of the GGL priorparameters (α_(o,j), β_(o,j))o∈

, 1≤j≤j_(max) from the corresponding detail coefficients (ζ _(o,j))o∈

, 1≤j≤j_(max).

On the one hand, ML estimators ({circumflex over (μ)}, {circumflex over(σ)}) are explicitly given by the empirical mean and standard deviation:

$\begin{matrix}{{{{\hat{\mu}}^{Re} = {\frac{1}{K_{j_{\max}}}{\sum\limits_{k = 1}^{K_{j_{\max}}}\;{{Re}\left( {\overset{\_}{\zeta}}_{a,k} \right)}}}},{{\hat{\sigma}}^{Re} = \sqrt{\frac{1}{K_{j_{\max}}}{\sum\limits_{k = 1}^{K_{j_{\max}}}\;\left( {{{Re}\left( {\overset{\_}{\zeta}}_{a,k} \right)} - {\hat{\mu}}^{Re}} \right)^{2}}}}}{{\hat{\mu}}^{Im} = {\frac{1}{K_{j_{\max}}}{\sum\limits_{k = 1}^{K_{j_{\max}}}\;{{Im}\left( {\overset{\_}{\zeta}}_{a,k} \right)}}}}{and}{{\hat{\sigma}}^{Im} = {\sqrt{\frac{1}{K_{j_{\max}}}{\sum\limits_{k = 1}^{K_{j_{\max}}}\;\left( {{{Re}\left( {\overset{\_}{\zeta}}_{a,k} \right)} - {\hat{\mu}}^{Im}} \right)^{2}}}.}}} & (25)\end{matrix}$

For each resolution level j and orientation o, {circumflex over(α)}_(o,j) ^(Re) and {circumflex over (β)}_(o,j) ^(Re) are estimatedfrom ζ _(o,j) as follows:

$\begin{matrix}\begin{matrix}{\left( {{\hat{\alpha}}_{o,j}^{Re},{\hat{\beta}}_{o,j}^{Re}} \right) = {\underset{{({\alpha,\beta})} \in {{\mathbb{R}}_{+} \times {\mathbb{R}}_{+}^{*}}}{\arg\mspace{11mu}\max}{f\left( {{{{Re}\left( {\overset{\_}{\zeta}}_{o,j} \right)};\alpha},\beta} \right)}}} \\{= {\underset{{({\alpha,\beta})} \in {{\mathbb{R}}_{+} \times {\mathbb{R}}_{+}^{*}}}{\arg\mspace{11mu}\max}{\sum\limits_{k = 1}^{K_{j}}\;{\log\mspace{11mu}{f\left( {{{{Re}\left( {\overset{\_}{\zeta}}_{o,j,k} \right)};\alpha},\beta} \right)}}}}} \\{= {\underset{{({\alpha,\beta})} \in {{\mathbb{R}}_{+} \times {\mathbb{R}}_{+}^{*}}}{\arg\mspace{11mu}\max}\left\{ {{\alpha{\sum\limits_{k = 1}^{K_{j}}\;{{{Re}\left( {\overset{\_}{\zeta}}_{o,j,k} \right)}}}} + {\frac{\beta}{2}{\sum\limits_{k = 1}^{K_{j}}\;{{Re}\left( {\overset{\_}{\zeta}}_{o,j,k} \right)}^{2}}} +} \right.}} \\\left. {\frac{K_{j}\alpha^{2}}{2\;\beta} - {\frac{K_{j}}{2}{\log\left( \frac{\beta}{2\pi} \right)}} + {K_{j}{\log\left( {{erfc}\left( \frac{\alpha}{\sqrt{2\;\beta}} \right)} \right)}}} \right\}\end{matrix} & (26)\end{matrix}$

The hyperparameters {circumflex over (α)}_(o,j) ^(Im) and {circumflexover (β)}_(o,j) ^(Im) are estimated in the same way, by replacing Re(·)by Im(·) in equation (26).

This two-dimensional minimization problem does not admit a closed formsolution. Hence, the ML parameters are computed using a numericaloptimization method, such as direct search method (eg, Rosenbrock'smethod, see [Bertsekas, 1995, chap 1, p 159-165]). Alternative solutionsbased on Monte Carlo methods or the Stein principle can also be thoughtof, at the expense of an increased computational burden.

When temporal regularization is considered, two additional parametershave to be determined: the temporal regularization parameter κ and theshape parameter p (see equations 21 and 22).

Indeed, in a subject-level analysis, only data recorded for one subjectis processed. For this reason, manually fixing the regularizationparameter κ and the shape parameter p may be suitable regarding to thereasonable amount of data, even if this makes the proposed method notcompletely automatic and remains somehow user-dependent. When agroup-level analysis is considered, the processed data becomes muchlarger, and since many subjects are considered, setting the temporalregularization parameter manually for all subjects at once is definitelysub-optimal due to the between-subject variability. This variability maybe significant either from an anatomical or functional viewpoint.Manually setting this parameter for each subject apart would bepossible, but cumbersome and it would make the proposed method moreuser-dependent regarding to the involved number of subjects in agroup-level analysis (typically 15 subjects involved).

Even in a group-level analysis, the estimation is performed for eachsubject apart. Moreover, a parameter value is preferably estimated foreach voxel of the 3D volume, which means that each voxel is processedindependently of the others. Nonetheless, for the sake of computationalefficiency, the same shape parameter p is considered for the wholevolume. Otherwise stated, instead of determining a unique temporalregularization parameter, κ, a respective parameter κ_(r) is determinedfor each voxel.

Again, a maximum likelihood estimation of said parameters is performed.Let ρ¹(r) . . . ρ^(Nr)(r) be the 3D images forming the 4D dataset, andlet define V_(r) as:V _(r)=(ρ²(r)−ρ¹(r),ρ³(r)−ρ²(r), . . . ,ρ^(2Nr)(r)−ρ^(Nr−1)(r))

Then, the joint estimation of the parameters p and κ_(r) in the ML senseis performed by minimizing the following criterion:L(p,κ)∝−log(κ_(r))+κ_(r)Σ_(t=1) ^(N) ^(r) ⁻¹ |V _(r)(t)|^(p)  (27)preferably under the constraint p≥1.

Obviously, a given voxel is supposed to more likely behave similarly toits neighbors, as well as for anatomical and functional perspective. Analternative method consists of estimating a parameter value perfunctional region of the brain (i.e. a group of neighboring voxelshaving similar functional roles). This estimation is based on an apriori classification step of the brain in order to identify thedifferent functional areas. This extension makes the estimation processmore robust since more samples are available at a given estimation step.Indeed, the temporal signals related to all the voxels belonging to thesame functional region are gathered in the same vector. If we denote byS, the size of a functional region t, the Maximum Likelihood estimationwill be made on a signal involving S, times more samples. The resultingestimation is therefore more robust to atypical observations (outliers),which may be occur in an fMRI session.

The image reconstruction method of [Chaâri et al. 2008] and [Chaâri etal. 2009] is based on the decomposition of the observed volume on awavelet basis. A further improvement, constituting another aspect of thepresent invention, consists in generalizing said method to the casewhere decomposition is performed on a frame.

Frames generalize the concepts of wavelets, curvelets, bandlets, etc. Aframe is a family of functions on which a signal (or image, which issimply a multi-dimensional signal) can be decomposed. In general, theresulting representation of the signal is redundant; a frame is called abasis when it is not redundant. Frames are generally used to helpcapturing more geometrical details of images under investigation. Forthis reason, using frames in pMRI reconstruction is helpful in order totrack reconstruction artefacts'.

More precisely, let us consider real-valued digital signals of length Las elements of the Euclidean space

^(L) endowed with the usual scalar product and norm denoted as

|

and ∥·∥, respectively.

Let K be an integer greater than or equal to L. A family of vectors(e_(k))_(1≤k≤K) in the finite-dimensional space

^(K) is a frame when there exists a constant μ in]0, +∞[ such that

$\begin{matrix}{{\forall{y \in {\mathbb{R}}^{L}}},{{\mu{y}^{2}} \leq {\sum\limits_{k = 1}^{K}\;{\left\langle y \middle| e_{k} \right\rangle }^{2}}}} & (28)\end{matrix}$

If the inequality (28) becomes an equality, is called a (e_(k))_(1≤k≤K)tight frame. The bounded linear frame analysis operator F and theadjoint synthesis frame operator F* are defined as

$\begin{matrix}{{F:\left. {\mathbb{R}}^{L}\rightarrow{{\mathbb{R}}^{K}:\left. y\mapsto\left( \left\langle y \middle| e_{k} \right\rangle \right)_{1 \leq k \leq K} \right.} \right.}{F^{*}:\left. {\mathbb{R}}^{K}\rightarrow{{\mathbb{R}}^{L}:\left. \left( \xi_{k} \right)_{1 \leq k \leq K}\mapsto{\sum\limits_{k = 1}^{K}\;{\xi_{k}e_{k}}} \right.} \right.}} & (29)\end{matrix}$

Note that F is injective whereas F* is surjective. When F¹=F*, then(e_(k))_(1≤k≤K) is an orthonormal basis. A simple example of a redundantframe is the union of orthonormal bases. In this case, the frame istight with μ=M and thus F*F=MI, where I is the identity operator.

In the following, two frame estimation problems of increasing complexitywill be addressed, the first one corresponding to a denoisingformulation and the second one to an extension involving a degradationoperator such as the sensitivity matrix S in parallel MRIreconstruction. The denoising problem is addressed first and itsapplication to pMRI matches the situation in which a noisy referenceimage, denoted as y in what follows, (e.g. such as the SENSEreconstruction) is available to estimate the hyperparameters of theframe representation. The second problem will be discussed afterpresenting Algorithm 4 and would correspond to estimating thehyperparameters of the frame decomposition directly from the reduced FOVimages in the pMRI context.

An observed signal y=F*x+n can be written according to its framerepresentation (FR) involving coefficients x∈

^(K) as follows:y=F*x+n  (30)where n is the error between the observed signal and its FR F*x. Thiserror is modeled by imposing that belongs to the closed convex setC _(δ) ={x∈

^(K) |N(y−F*x)≤δ}

where δ∈[0, ∞[ is some error bound and N can be any norm on

^(L).

In signal/image recovery problems, n is an additive noise that corruptsthe measured data. The following developments will be focused on thecase of a bounded observation error modeled by uniform noise. Byadopting a probabilistic approach, y and x are assumed to berealizations of random vectors Y and X. In this context, the goal is tocharacterize the probability distribution of XIY by considering someparametric probabilistic model and by estimating the associatedhyperparameters. The estimation problem is much more difficult than inthe case of the decomposition on a non-redundant basis, e.g. a waveletbasis, because F. is not bijective.

In a Bayesian framework, it is necessary to define prior distributionsfor the frame coefficients. For instance, this prior may be chosen so asto promote the sparsity of the representation. In the following f(x,θ),denotes the probability density function (pdf) of the frame coefficientsthat depends on an unknown hyperparameter vector θ and f(θ), is the apriori pdf for the hyperparameter vector θ. In compliance with theconstraint C_(δ), n is assumed to be uniformly distributed on the ballB_(0,δ)={α∈

^(L)|N(α)≤δ}, and ƒ(y|x) is the uniform pdf on the closed convex ball:B _(F*x,δ) ={y∈

^(L) |N(y−F*x)≤δ}

The hyperparameter vector θ is also considered a realization of therandom vector Θ, and the conditional pdf of (X,Θ) given Y can be writtenas:ƒ(x,θ|y)∝ƒ(y|x)ƒ(x|θ)ƒ(θ)  (31)

For the sake of simplicity, it will be assumed that frame coefficientsare a priori independent with marginal GG (Generalized Gaussian)distributions even if, as discussed above, a GGL distribution is bettersuited for applications to pMRI This leads to the following framecoefficient prior:

$\begin{matrix}{{f\left( {\left. x_{k} \middle| \alpha_{k} \right.,\beta_{k}} \right)} = {\frac{\beta_{k}}{2\alpha_{k}{\Gamma\left( \frac{1}{\beta_{k}} \right)}}{\exp\left( {- \frac{{x_{k}}^{\beta_{k}}}{\alpha_{k}^{\beta_{k}}}} \right)}}} & (32)\end{matrix}$

where α_(k)>0, β_(k)>0 (with k∈{1, . . . , K}) are the scale and shapeparameters associated with x_(k), which is the component of the framecoefficient vector and Γ(·) is the Gamma function.

By introducing it γ_(k)=α_(k) ^(β) ^(k) , the frame prior can berewritten as:

$\begin{matrix}{{f\left( {\left. x_{k} \middle| \gamma_{k} \right.,\beta_{k}} \right)} = {\frac{\beta_{k}}{2\gamma_{k}^{1/\beta_{k}}{\Gamma\left( \frac{1}{\beta_{k}} \right)}}{\exp\left( {- \frac{{x_{k}}^{\beta_{k}}}{\gamma_{k}}} \right)}}} & (33)\end{matrix}$

The distribution of a frame coefficient generally differs from onecoefficient to another. However, some frame coefficients can have verysimilar distributions, which can be defined by the same hyperparameters.As a consequence, it is proposed to split the frame coefficients into Gdifferent groups. The g^(th) group will be parameterized by a uniquehyperparameter vector denoted as {dot over (θ)}_(g)=(β_(g), γ_(g)). Inthis case, the frame prior can be expressed as

$\begin{matrix}{{f\left( x \middle| \theta \right)} = {\prod\limits_{g = 1}^{G}\;\left\lbrack {\left( \frac{\beta_{g}}{2\gamma_{g}^{1/\beta_{s}}{\Gamma\left( \frac{1}{\beta_{g}} \right)}} \right)^{n_{g}}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)}} \right\rbrack}} & (34)\end{matrix}$where the summation covers the index set S_(g) of the elements of theg^(th) group containing n_(g) elements and θ_(g)=(θ₁ . . . θ_(G)). E.g.,each group can correspond to a given wavelet subband. A coarserclassification may be mad when using multiscale frame representations byconsidering that all the frame coefficients at a given resolution levelbelong to a same group.

The hierarchical Bayesian model for the frame decomposition is completedby the following improper hyperprior:

$\begin{matrix}{{f(\theta)} = {{\prod\limits_{g = 1}^{G}\;{f\left( \theta_{g} \right)}} = {{\prod\limits_{g = 1}^{G}\;\left\lbrack {{f\left( \gamma_{g} \right)}{f\left( \beta_{g} \right)}} \right\rbrack} \propto {\prod\limits_{g = 1}^{G}\;\left\lbrack {{\frac{1}{\gamma_{g}}1_{\mathbb{R}}} + {\left( \gamma_{g} \right)1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)}} \right\rbrack}}}} & (35)\end{matrix}$where 1_(A)(ξ) is the function defined on A⊂

by 1_(A)(ξ)=1 if ξ∈A and 1_(A)(ξ)=0 otherwise.

The motivations for using this kind of prior are the following

-   -   the interval [0,3] covers all possible values of β_(g)        encountered in practical applications, and there is no        additional information about the parameter.    -   The prior for the parameter γ_(g) is a Jeffrey's distribution        that reflects the absence of knowledge about this parameter.

The resulting posterior distribution is therefore given by:

$\begin{matrix}{{f\left( {x,\left. 0 \middle| y \right.} \right)} \propto {1_{C_{\delta}}(x){\quad{\underset{\;{g = 1}}{\overset{G}{\prod}}{\quad\;\left\lbrack {\left( \frac{\beta_{g}}{2\gamma_{g}^{1/\beta_{g}}{\Gamma\left( \frac{1}{\beta_{g}} \right)}} \right)^{n_{g}}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{0}}\;{x_{k}}^{\beta_{o}}}} \right)}\left( {\frac{1}{\gamma_{g}}1_{{\mathbb{R}} +}\left( \gamma_{g} \right)1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)} \right)} \right\rbrack}}}}} & (36)\end{matrix}$

The Bayesian estimators [e.g., the maximum a posteriori (MAP) or minimummean square error (MMSE) estimators also known as the posterior meanestimator] associated with the posterior distribution (35) have nosimple closed-form expression.

In what follows, a stochastic procedure is proposed to compute the MMSEestimator, which relies on hybrid Markov Chain Monte Carlo (MCMC)algorithms. The idea is to use a suitable algorithm to generate samplesdistributed according to the posterior distribution (35). Afterconvergence, the generated samples are used to compute the MMSEestimates of the unknown model parameter and hyperparameter vectors xand θ, respectively. In the pMRI context, x is the sought framerepresentation of the reference image.

When the considered frame is the union of orthonormal M bases and N(·)is the Euclidean norm, the well-known Gibbs sampler (GS) can be used[Geman, 1984], which iteratively generates samples distributed accordingto conditional distributions associated with the target distribution.More precisely, the basic GS iteratively generates samples distributedaccording to f(x|θ,y) and f(θ|x, y) to simulate realizations of the fullposterior. f(x, θ|y).

Straightforward calculations yield the following conditionaldistribution:

$\begin{matrix}{{f\left( {\left. x \middle| \theta \right.,y} \right)} \propto {1_{C_{\delta}}(x){\prod\limits_{g = 1}^{G}\;{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)}}}} & (37)\end{matrix}$

This conditional distribution is a product of GG distributions truncatedon C. Actually, sampling according to this truncated distribution is notalways easy to perform since the adjoint frame operator F_is usually oflarge dimension. However, two alternative sampling strategies aredetailed in what follows.

Naive sampling proceeds by sampling according to independent GGdistributions

$\prod\limits_{g = 1}^{G}\;{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)}$and then accepting the proposed candidate x only if N(y−F*x)≤δ. Thismethod can be used for any frame decomposition and any norm. However, itcan be quite inefficient because of a very low acceptance ratio,especially when δ takes small values.

The Gibbs sampler is designed to sample more efficiently from theconditional distribution in Eq. (36) when the considered frame is theunion of M orthonormal bases and N(·) is the Euclidean norm. In thiscase, the analysis frame operator and the corresponding adjoint can bewritten as F=[F₁ . . . F_(M)]^(T) and F*=[F*₁ . . . F*_(M)]^(T)respectively, where ∀m∈{1, . . . , M}, F_(m) is the decompositionoperator onto the m^(th) orthonormal basis such as F*_(m)F_(m)=F_(m)F*_(m)=Id. In what follows, every x∈

^(K) with K=ML is decomposed as x=x₁ ^(T), . . . , x_(M) ^(t)]^(T) wherex_(m)∈

^(L) ∀m∈{1, . . . , M}.

The GS for the generation of frame coefficients draws vectors accordingto the conditional distribution ƒ(x_(n)|x_(−n), y, θ) under theconstraint N(y−F*x)≤δ where x_(−n) is the reduced size vector ofdimension

^(K−L) built from x by removing the n^(th) vector x_(n). If N(·) is theEuclidean norm, then ∀n∈{1, . . . , M}:

$\begin{matrix}{\left. {{N\left( {y - {\sum\limits_{m = 1}^{M}\;{F_{m}^{*}x_{m}}}} \right)} \leq \delta}\Leftrightarrow{{{F_{n}^{*}\left( {{F_{n}y} - {\sum\limits_{m = 1}^{M}\;{F_{n}F_{n}^{*}x_{m}}}} \right)}} \leq \delta}\Leftrightarrow{{{{F_{n}y} - {\sum\limits_{m \neq n}\;{F_{n}F_{m}^{*}x_{m}}} - x_{n}}} \leq {\delta\left( {{{since}\mspace{14mu}{\forall{z \in {\mathbb{R}}^{L}}}},{{{F_{n}^{*}z}} = {z}}} \right)}}\Leftrightarrow{{N\left( {x_{n} - c_{n}} \right)} \leq \delta} \right.,\mspace{20mu}{{{where}\mspace{14mu} c_{n}} = {{F_{n}y} - {\sum\limits_{m \neq n}\;{F_{n}F_{m}^{*}x_{m}}}}}} & (38)\end{matrix}$

To sample each x_(n), it is proposed to use an MH step whose proposaldistribution is supported on the ball B_(Cn,δ) defined by:B _(Cn,δ)={α∈

^(L) |N(α−c _(n))≤δ}  (39)

Random generation from a pdf q_(δ) defined on B_(Cn,δ) will now bediscussed.

First of all, it will be considered how to sample vectors in the unitl_(p) ball (p∈]0,+∞]) of

^(L). In the special case p=+∞, this can be easily performed by samplingindependently along each space coordinate according to a distribution onthe interval [−1, 1]. The parameter “p” considered here should not beconfused with the shape parameter used in temporal regularization, seeequations 21 and 22.

Thus, this section focuses on the more difficult problem associated witha finite value of p. In the following, ∥·∥_(p) denotes the l_(p) norm.

Let A=[A₁, . . . , A_(L′)]^(T) be the random vector of i.i.d. componentswhich have the following GG(P^(1/p), p) pdf:

${\forall{a \in {\mathbb{R}}}},{{f(a)} = {\frac{p^{1 - {1/p}}}{2{\Gamma\left( {1/p} \right)}}{\exp\left( {- \frac{{a}^{p}}{p}} \right)}}}$

Let U=[U₁, . . . , U_(L′)]^(T)=A/∥A∥_(p). Then, it can be shown ([Guptaet al., 1997]) that the random vector U is uniformly distributed on thesurface of the l_(p) unit sphere of

^(L′) and the joint pdf of U₁, . . . , U_(L′−1) is

${f\left( {u_{1},\ldots\mspace{14mu},u_{L^{\prime} - 1}} \right)} = {\frac{p^{L^{\prime} - 1}{\Gamma\left( {L^{\prime}/p} \right)}}{2^{L^{\prime} - 1}\left( {\Gamma\left( {1/p} \right)} \right)^{L^{\prime}}}\left( {1 - {\sum\limits_{k = 1}^{L^{\prime} - 1}\;{u_{k}}^{p}}} \right)^{{({1 - p})}/p}1_{D_{p,L^{\prime}}}\left( {u_{1},\ldots\mspace{14mu},u_{L^{\prime} - 1}} \right)}$$\mspace{79mu}{{{where}\mspace{20mu} D_{p,L^{\prime}}} = {\left\{ {\left( {u_{1},\ldots\mspace{14mu},u_{L^{\prime} - 1}} \right) \in {{\mathbb{R}}^{L^{\prime} - 1}{\sum\limits_{k = 1}^{L^{\prime} - 1}\; }u_{k}}} \middle| {}_{p}{< 1} \right\}\mspace{11mu}.}}$

The uniform distribution on the unit l_(p) sphere of

^(L′) will be denoted by U(L′,p). Therefore U=[U₁, . . . ,U_(L′)]^(T)˜U(L′,p).

For every ∈{1, . . . , L′−1}, it can be shown ([Gupta et al., 1997])that the pdf of V=[U₁, . . . , U_(L)]^(T) is given by:

$\begin{matrix}{{q_{1}\left( {u_{1},\ldots\mspace{14mu},u_{L}} \right)} = {\frac{p^{L}{\Gamma\left( {L^{\prime}/p} \right)}\left( {1 - {\sum\limits_{k = 1}^{L}\;{u_{k}}^{p}}} \right)^{{{({L^{\prime} - L})}/p} - 1}}{2^{L}\left( {\Gamma\left( {1/p} \right)} \right)^{L}{\Gamma\left( {\left( {L^{\prime} - L} \right)/p} \right)}}1_{D_{p,{L + 1}}}\left( {u_{1},\ldots\mspace{14mu},u_{L}} \right)}} & (40)\end{matrix}$

In particular, if p∈

* and L′=L+p, equation (40) provides the uniform distribution on theunit l_(p) sphere of

^(L′). Sampling from any distribution q_(η) on the ball of radius δ≥η>0is straightforwardly deduced by scaling V.

Having a closed form expression of this pdf is important to be able tocalculate the acceptance ratio of the MH move. To take into account thevalue of x_(n) ^((i−1)) obtained at the previous iteration (i−1), it mayhowever be preferable to choose a proposal distribution supported on arestricted ball of radius η∈]0, δ[containing x_(n) ^((i−1)). Thisstrategy, similar to the random walk MH, results in a better explorationof regions associated with large values of the conditional distributionƒ(x|θ, y). More precisely, it is proposed to choose a proposaldistribution defined on B_({circumflex over (x)}) _(n) _((i−1)) , ηwhere{circumflex over (x)} _(n) ^((i−1)) =P(x _(n) ^((i−1)) −c _(n))+c_(n)  (41)and P is the projection onto the ball B_(0,δ−η) defined as:

$\begin{matrix}{{\forall{a \in {\mathbb{R}}^{L}}},\mspace{14mu}{{P(a)} = \left\{ \begin{matrix}a & {{{if}\mspace{14mu}{N(a)}} \leq {\delta - \eta}} \\{\frac{\delta - \eta}{N(a)}a} & {{otherwise}.}\end{matrix} \right.}} & (42)\end{matrix}$

This choice of the center of the ball guarantees thatB_({circumflex over (x)}) _(n) _((i−1)) _(,η)⊂B_(c) _(n) _(,δ).

Moreover, any point of B_(c) _(n) _(,δ) can be reached after consecutivedraws in B_({circumflex over (x)}) _(n) _((i−1)) _(,η). The radius η hasto be adjusted to ensure a good exploration of B_(c) _(n) _(,δ). Inpractice, it may also be interesting to fix a small enough value of η(compared with δ) so as to improve the acceptance ratio.

Instead of sampling the hyper-parameter vector θ according to ƒ(θ|x, y)it is proposed to iteratively sample according to ƒ(γ_(g)|β_(g), x, y)and ƒ(β_(g)|γ_(g), x, y). Straightforward calculations yield thefollowing results:

$\begin{matrix}{{f\left( {\left. \gamma_{g} \middle| \beta_{g} \right.,x,y} \right)} \propto {{\gamma_{g}^{{- \frac{n_{g}}{\beta_{g}}} - 1}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)}1_{\mathbb{R}}} + \left( \gamma_{g} \right)}} & (43) \\{{f\left( {\left. \beta_{g} \middle| \gamma_{g} \right.,x,y} \right)} \propto {\frac{\beta_{g}^{n_{g}}1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)}{{\gamma_{g}^{n_{g}/\beta_{g}}\left\lbrack {\Gamma\left( {1/\beta_{g}} \right)} \right\rbrack}^{n_{g}}}{{\exp\left( {\sum\limits_{k \in S_{g}}\;\frac{- {x_{k}}^{\beta_{g}}}{\gamma_{g}}} \right)}.}}} & (44)\end{matrix}$

Consequently, ƒ(γ_(g)|β_(g), x, y) is the pdf of the inverse gammadistribution

${\mathcal{I}\mathcal{G}}\left( {\frac{n_{g}}{\beta_{g}},{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)$that is easy to sample. Conversely, it is more difficult to sampleaccording to the truncated pdf ƒ(β_(g)|γ_(g), x, y). This is achieved byusing an MH move whose proposal q(β_(g)|β_(g) ^((i−1))) is a Gaussiandistribution truncated on the interval [0, 3] with standard deviationσβ_(a)=0.05 [Dobigeon et al, 2007], [Robert, 1995]. The mode of thisdistribution is the value of the parameter β_(g) ^((i−1)) at theprevious iteration (i−1). The resulting method is the hybrid GSsummarized in Algorithm 4.

If one wants to integrate altogether the reconstruction andhyperparameter estimation steps in the same Bayesian framework,observation model in Eq. 30 can be extended to y=SF*x+n, where S is thesensitivity matrix defined in Eq. (3). The inference has then to bereconducted accordingly.

Algorithm 4 1: Initialize with some θ⁽⁰⁾ = (θ_(g) ⁽⁰⁾)1 ≤ g ≤ G = (γ_(g)⁽⁰⁾, β_(g) ⁽⁰⁾)_(1≤g≤G) and x⁽⁰⁾ ∈ C_(δ), and set i = 1. 2: repeat 3:Sampling x: 4: for n = 1 to M do 5: Compute c_(n) ^((i)) = F_(n)(y −Σ_(m<n) F_(m) * x_(m) ^((i)) − Σ_(m>n) F_(m) * x_(m) ^((i−1))) and{circumflex over (x)}_(n) ^((i−1)) = P(x_(n) ^((i−1)) − c_(n) ^((i))) +c_(n) ^((i)). 6: Simulate x_(n) ^((i)) as follows: Generate {tilde over(x)}_(n) ^((i)) ~ q_(η)(x_(n) − {circumflex over (x)}_(n) ^((i−1)))where q_(η) is defined on β_(o,η) Compute the ratio${r\left( {{\overset{\sim}{x}}_{n}^{(i)},x_{n}^{({i - 1})}} \right)} = \frac{{f\left( {{{\overset{\sim}{x}}_{n}^{(i)}❘\theta^{({i - 1})}},{{\left( x_{m}^{(i)} \right)m} < n},{{\left( x_{m}^{({i - 1})} \right)m} > n},y} \right)}{q_{n}\left( {x_{n}^{({i - 1})} - {P\left( {{\overset{\sim}{x}}_{n}^{(i)} - c_{n}^{(i)}} \right)} - c_{n}^{(i)}} \right)}}{{f\left( {{x_{n}^{({i - 1})}❘\theta^{({i - 1})}},{{\left( x_{m}^{(i)} \right)m} < n},{{\left( x_{m}^{({i - 1})} \right)m} > n},y} \right)}{q_{n}\left( {{\overset{\sim}{x}}_{n}^{(i)} - {\hat{x}}_{n}^{({i - 1})}} \right)}}$and accept the proposed candidate with the probability min{1, r({tildeover (x)}_(n) ^((i)), x_(n) ^((i−1)))} 7: end for 8: Sampling θ: 9: forg = 1 to G do 10:${Generate}\mspace{14mu}{\left. \gamma_{g}^{(i)} \right.\sim I}\;{G\left( {\frac{n_{g}}{\beta_{g}^{({i - 1})}},{\sum\limits_{k \in S_{g}}{❘{x_{k}^{(i)}❘\beta_{g}^{({i - 1})}}}}} \right)}$11: Simulate β_(g) ^((i)) as follows: Generate {tilde over (β)}_(g)^((i)) ~ q(β_(g)|β_(g) ^((i−1))). Compute the ratio${r\left( {{\overset{\sim}{\beta}}_{g}^{(i)},\beta_{g}^{({i - 1})}} \right)} = \frac{{f\left( {{{\overset{\sim}{\beta}}_{g}^{(i)}❘\gamma_{g}^{(i)}},x^{(i)},y} \right)}{q\left( {\beta_{g}^{({i - 1})}❘{\overset{\sim}{\beta}}_{g}^{(i)}} \right)}}{{f\left( {{\beta_{g}^{({i - 1})}❘\gamma_{g}^{(i)}},x^{(i)},y} \right)}{q\left( {{\overset{\sim}{\beta}}_{g}^{(i)}❘\beta_{g}^{({i - 1})}} \right)}}$and accept the proposed candidate with the probability min{1, r({tildeover (β)}_(g) ^((i)), β_(g) ^((i−1)))}. 12: end for 13: Set i ← t + 1.14: Until Convergence

Although this algorithm is intuitive and simple to implement, it must bepointed out that it was derived under the restrictive assumption thatthe considered frame is the union of M orthonormal bases. When theseassumptions do not hold, another algorithm, discussed below, allowssampling frame coefficients and the related hyperparameters byexploiting algebraic properties of frames. Indeed, as a directgeneration of samples according to ƒ(x|θ, y) is generally impossible, analternative method is proposed which replaces the Gibbs move by an MHmove. This MH move aims at sampling globally a candidate x according toa proposal distribution. This candidate is accepted or rejected with thestandard MH acceptance ratio. The efficiency of the MH move stronglydepends on the choice of the proposal distribution for x.

We denote as x^((i)) be the i^(th) accepted sample of the algorithm andq(x|x^((i−1))) the proposal that is used to generate a candidate atiteration i. The main difficulty for choosing q(x|x^((i−1))) stems fromthe fact that it must guarantee that x∈C_(δ) while yielding a tractableexpression of q(x^((i−1))|x)/q(x|x^((i−1))).

For this reason, it is proposed to exploit the algebraic properties offrame representations. More precisely, any frame coefficient vector canbe decomposed as x=x_(H)+x_(⊥), where x_(H) and x_(H⊥) are realizationsof random vectors taking their values in H=Ran(F) and H^(⊥),=[Ran(F)]^(⊥)=Null(F*), respectively. It is recalled that the range of Fis Ran(F)={x∈

^(K)|∃y∈

^(L), Fy=x} and the null space of F* isNull(F*)={x∈

^(K) |F*x=0}

The proposal distribution used here allows generating samples x_(H)∈Hand x_(H⊥), ∈H^(⊥). More precisely, the following separable form of theproposal pdf will be considered:q(x|x ^((i)))=q(x _(H) |x _(H) ^((i−1)))q(x _(H) _(⊥) |x _(H) _(⊥)^((i−1)))  (45).)where x_(H) ^((i−1))∈H, x_(H) _(⊥) ^((i−1))∈H^(⊥) and x^((i−1)=x) _(H)^((i−1))+x_(H) _(⊥) ^((i−1)). In other words, independent sampling ofx_(H) and x_(H⊥) will be performed.

If the decomposition x=x_(H)+x_(H⊥) is considered, sampling x in C_(δ)is equivalent to sampling λ∈C _(δ), where C _(δ)={λ∈

^(L)|N(y−F*Fλ)≤δ}. Indeed, x_(H)=Fλ where λ∈

^(L) and, since x_(H) _(⊥) ∈Null(F*), F*x=F*Fλ. Sampling λ in C _(δ) canbe easily achieved, e.g., by generating u from a distribution on theball B_(y,δ) and by taking λ=(F*F)⁻¹ u.

To make the sampling of x_(H) at iteration i more efficient, taking intoaccount the sampled value at the previous iteration x_(H)^((i−1))=Fλ^((i−1))=F(F*F)⁻¹u^((i−1)) may be interesting.

Similarly to random walk generation techniques, u is generated in B_(û)_((i−1)) _(,η) where η∈]0, δ[ and û^((i−1))=P(u^((i−1))−y)+y. Thisallows drawing a vector u such that x_(H)=F(F*F)⁻¹u∈C_(δ) andN(u−u^((i−1)))≤2η.

The generation of u can then be performed according to equation (40)provided that N(·) is an l_(p) norm with p∈[1,+∞].

Once that x_(H)=Fλ∈H∩C_(δ): (which ensures that x is in C_(δ)), has beensimulated, x_(H⊥) has to be sampled as an element of H^(⊥). Sincey=n=F*x+n=F*x_(H)+n, there is no information in y about x_(H⊥). As aconsequence, it is proposed to sample x_(H) by drawing z according tothe Gaussian distribution

(x^((i−1)), σ_(x) ²I) and by projecting z onto H^(⊥), i.e.,x _(H) _(⊥) =Π_(H) _(⊥) z

Π_(H) _(⊥) =I−F(F*F)⁻¹F* is the orthogonal projection operator ontoH^(⊥).

It can then be shown that the expression of the proposal pdf

$\begin{matrix}{\frac{q\left( x^{({i - 1})} \middle| x \right)}{q\left( x \middle| x^{({i - 1})} \right)} = {\frac{q_{\eta}\left( {u^{({i - 1})} - {P\left( {u - y} \right)} - y} \right)}{q_{\eta}\left( {u - {\hat{u}}^{({i - 1})}} \right)}.}} & (46)\end{matrix}$

This expression remains valid in the degenerate case when K=L (yieldingx_(H⊥)=0). Finally, it is important to note that, if q_(η) is theuniform distribution on the ball B_(0, η), the above ratio reduces to 1,which simplifies the computation of the MH acceptance ratio. The finalalgorithm is summarized in Algorithm 5. It is noteworthy that thesampling of the hyper-parameter vector is performed as for the hybrid GSof Algorithm 4.

Algorithm 5 1: Initialize with some θ⁽⁰⁾ = (θ_(g) ⁽⁰⁾)_(1≤g≤G) = (γ_(g)⁽⁰⁾, β_(g) ⁽⁰⁾)_(1≤g≤G) and u⁽⁰⁾ ∈ B_(y,δ). Set x⁽⁰⁾ = F(F * F)⁻¹u⁽⁰⁾and i = 1. 2: Repeat 3: Sampling x: Compute û^((i−1)) = P(u^((i−1)) −y) + y. Generate ũ^((i)) ~ q_(n)(u − û^((i−1))) where q_(n) is definedon β_(o,η). Compute {tilde over (x)}_(H) ^((i)) = F(F * F)⁻¹ũ^((i)).Generate z^((i)) ~ N(x^((i−1)),σ_(x) ²I). Compute {tilde over (x)}_(H)_(⊥) ^((i)) = Π_(H) ^(⊥)z^((i)) and {tilde over (x)}^((i)) = {tilde over(x)}_(H) ^((i)) + {tilde over (x)}_(H) ^(⊥(i)).${{Compute}\mspace{14mu}{the}\mspace{14mu}{ratio}\mspace{14mu} r\left( {{\overset{\sim}{x}}^{(i)},x^{({i - 1})}} \right)} = \frac{{f\left( {{{\overset{\sim}{x}}^{(i)}❘\theta^{({i - 1})}},y} \right)}{q_{n}\left( {u^{({i - 1})} - {P\left( {{\overset{\sim}{u}}^{(i)} - y} \right)} - y} \right)}}{{f\left( {{x^{({i - 1})}❘\theta^{({i - 1})}},y} \right)}{q_{n}\left( {{\overset{\sim}{u}}^{(i)} - {\overset{\sim}{u}}^{({i - 1})}} \right)}}$and accept the proposed candidates ũ^((i)) and {tilde over (x)}^((i))with probability min{1, r({tilde over (x)}^((i)), x^((i−1)))}. 4:Sampling θ: 5: for g = 1 to G do 6:${Generate}\mspace{14mu}{\left. \gamma_{g}^{(i)} \right.\sim I}\;{{G\left( {\frac{n_{g}}{\beta_{g}^{({i - 1})}},{\sum\limits_{k \in S_{g}}{❘{x_{k}^{(i)}❘\beta_{g}^{({i - 1})}}}}} \right)}.}$7: Simulate β_(g) ^((i)) as follows: Generate β_(g) ^((i)) ~q(β_(g)|β_(g) ^((i−1))). Compute the ratio${r\left( {{\overset{\sim}{\beta}}_{g}^{(i)},\beta_{g}^{({i - 1})}} \right)} = \frac{{f\left( {{{\overset{\sim}{\beta}}_{g}^{(i)}❘\gamma_{g}^{(i)}},x^{(i)},y} \right)}{q\left( {\beta_{g}^{({i - 1})}❘{\overset{\sim}{\beta}}_{g}^{(i)}} \right)}}{{f\left( {{\beta_{g}^{({i - 1})}❘\gamma_{g}^{(i)}},x^{(i)},y} \right)}{q\left( {{\overset{\sim}{\beta}}_{g}^{(i)}❘\beta_{g}^{({i - 1})}} \right)}}$and accept the proposed candidate with the probability min{1, r({tildeover (β)}_(g) ^((i)), β_(g) ^((i−1)))}. 8: end for 9: Set i ← i + 1. 10:Until Convergence

The hierarchical Bayesian model of equations (31-35) can be extended toinclude an additional term in the prior depending on the Total Variation(TV) of the image to be reconstructed. Like redundant framerepresentations, using TV priors leads also to a hyper-parameterestimation problem.

Assuming an exponential shape, the new prior can be expressed as:

$\begin{matrix}{{f\left( x \middle| \theta \right)} = {\frac{1}{Z(\theta)}{\exp\left( {{- \kappa}{{F^{*}x}}{TV}} \right)}{\prod\limits_{g = 1}^{G}\;\left\lbrack {\left( \frac{1}{\gamma_{g}^{1/\beta_{g}}} \right)^{n_{g}}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{s}}}} \right)}} \right\rbrack}}} & (47)\end{matrix}$

where θ=(θ₁, . . . , θ_(G), κ) is the new hyper-parameter vector withκ>0, ∥·∥_(TV) is the TV semi-norm and Z(θ) is a normalization constant.The new hierarchical Bayesian model for the frame decomposition iscompleted by the following improper hyperprior:

$\begin{matrix}{{{f(\theta)} = {{{Z(\theta)}{\varphi(\kappa)}{\prod\limits_{g = 1}^{G}\;{\varphi\left( \theta_{g} \right)}}} = {{{Z(\theta)}{\varphi(\kappa)}{\prod\limits_{g = 1}^{G}\;\left\lbrack {{\varphi\left( \gamma_{g} \right)}{\varphi\left( \beta_{g} \right)}} \right\rbrack}} \propto {{Z(\theta)}1_{\lbrack{0,\kappa_{\max}}\rbrack}(\kappa){\prod\limits_{g = 1}^{G}\;\left\lbrack {{\frac{1}{\gamma_{g}}1_{\mathbb{R}}} + {\left( \gamma_{g} \right)1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)}} \right\rbrack}}}}},} & (48)\end{matrix}$

where κ_(max) is a positive real to be fixed (in a preferred embodiment,κ_(max)=10).

It can be shown that Z(θ) is uniformly bounded with respect to γ_(g) andtherefore that the hyperprior f(θ) in Eq. (47) has a stable asymptoticbehaviour when γ_(g)→+∞.

The resulting new posterior distribution is therefore given by:

$\begin{matrix}{{f\left( {x,\left. \theta \middle| y \right.} \right)} = {1_{C^{*}}(x){\prod\limits_{g = 1}^{G}\;{\left\lbrack {\left( \frac{1}{\gamma_{g}^{1/\beta_{g}}} \right)^{n_{g}}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{s}}}} \right)}\left( {{\frac{1}{\gamma_{g}}1_{\mathbb{R}}} + {\left( \gamma_{g} \right)1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)}} \right)} \right\rbrack \times {\exp\left( {{- \kappa}{{F^{*}x}}{TV}} \right)}1_{\lbrack{0,\kappa_{\max}}\rbrack}{(\kappa).}}}}} & (49)\end{matrix}$

The Bayesian estimators associated with the posterior distribution inEq. (48) still have no simple closed-form expression. For this reason,the same sampling strategy discussed above will be applied, and theframe coefficients will be sampled as in Algorithm 5. However, for thehyper-parameters vector, straightforward calculations show that theposterior distribution for the hyper-parameters γ_(g), β_(g) and κ willbe expressed as:

$\begin{matrix}{\mspace{79mu}{{f\left( {\left. \gamma_{g} \middle| \beta_{g} \right.,\kappa,x,y} \right)} \propto {{\gamma_{g}^{{- \frac{n_{g}}{\beta_{g}}} - 1}{\exp\left( {{- \frac{1}{\gamma_{g}}}{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)}1_{\mathbb{R}}} + \left( \gamma_{g} \right)}}} & (50) \\{\mspace{79mu}{{{f\left( {\left. \beta_{g} \middle| \gamma_{g} \right.,\kappa,x,y} \right)} \propto {{\exp\left( {\sum\limits_{k \in S_{g}}\;{{- \frac{1}{\gamma_{g}}}{x_{k}}^{\beta_{g}}}} \right)}1_{\lbrack{0,3}\rbrack}\left( \beta_{g} \right)}}\mspace{20mu}{and}}} & (51) \\{{f\left( {\left. \kappa \middle| \gamma_{1} \right.,{\ldots\mspace{14mu}\gamma_{G}},\beta_{1},\ldots\mspace{14mu},\beta_{G},x,y} \right)} \propto {{\exp\left( {{- \kappa}{{F^{*}x}}{TV}} \right)}1_{\lbrack{0,\kappa_{\max}}\rbrack}(\kappa)}} & (52)\end{matrix}$respectively.

Consequently, ƒ(γ_(g)|β_(g), κ, x, y) is the pdf of the inverse gamma

${\mathcal{I}\mathcal{G}}\left( {\frac{n_{g}}{\beta_{g}},{\sum\limits_{k \in S_{g}}\;{x_{k}}^{\beta_{g}}}} \right)$Sampling γ_(g) will therefore be performed exactly as in Algorithm 4.Conversely, it is more difficult to sample according to ƒ(β_(g)|γ_(g),κ, x, y) and ƒ(κ|γ₁, . . . , γ_(G), β₁, . . . , β_(G), x, y). This taskis achieved by using two MH moves whose proposal distributionsq(β_(g)|β_(g) ^((i−1))) and q(κ|κ^((i−1))) are Gaussian distributionstruncated on the intervals [0, 3] and [0, κ_(max)] with standarddeviations σ_(βg)=0.05 and σ_(κ)=0.01, respectively. These standarddeviation values have been fixed based on our practical observations.The resulting method to sample according to the posterior distributionin Eq. (49) is summarized in Algorithm 6.

Algorithm 6 1: Initialize with some θ⁽⁰⁾ = ((θ_(g) ⁽⁰⁾)_(1≤g≤G), κ⁽⁰⁾) =((γ_(g) ⁽⁰⁾, β_(g) ⁽⁰⁾)_(1≤g≤G), κ⁽⁰⁾) and u⁽⁰⁾ ∈ B_(y,δ). Set x⁽⁰⁾ =F(F * F)⁻¹u⁰⁾ and i = 1. 2: repeat 3: Sampling x: Compute û^((i−1)) =P(u^((i−1)) − y) + y. Generate ũ^((i)) ~ q_(n)(u − û^((i−1))) whereq_(n) is defined on β_(o,η). Compute {tilde over (x)}_(H) ^((i)) =F(F*F)⁻¹ũ^((i)). Generate z^((i)) ~ N(x^((i−1)),σ_(x) ²I). Compute{tilde over (x)}_(H) ^(⊥(i)) = Π_(H) ^(⊥)z^((i)) and {tilde over(x)}^((i)) = {tilde over (x)}_(H) ^((i)) + {tilde over (x)}_(H) ^(⊥(i)).${{Compute}\mspace{14mu}{the}\mspace{14mu}{ratio}\mspace{14mu}{r\left( {{\overset{\sim}{x}}^{(i)},x^{({i - 1})}} \right)}} = \frac{{f\left( {{{\overset{\sim}{x}}^{(i)}❘\theta^{({i - 1})}},y} \right)}{q_{n}\left( {u^{({i - 1})} - {P\left( {{\overset{\sim}{u}}^{(i)} - y} \right)} - y} \right)}}{{f\left( {{x^{({i - 1})}❘\theta^{({i - 1})}},y} \right)}{q_{n}\left( {{\overset{\sim}{u}}^{(i)} - {\hat{u}}^{({i - 1})}} \right)}}$and accept the proposed candidates ũ^((i)) and {tilde over (x)}^((i))with probability min{1, r({tilde over (x)}^((i)), x^((i−1)))}. 4:Sampling θ: 5: for g − 1 to G do 6:${Generate}\mspace{14mu}{\left. \gamma_{g}^{(i)} \right.\sim I}\;{{G\left( {\frac{n_{g}}{\beta_{g}^{({i - 1})}},{\sum\limits_{k \in S_{g}}{❘{x_{k}^{(i)}❘\beta_{g}^{({i - 1})}}}}} \right)}.}$7: Simulate β_(g) ^((i)) as follows: Generate {tilde over (β)}_(g)^((i)) ~ q(β_(g)|β_(g) ^((i+1))).${{Compute}\mspace{14mu}{the}\mspace{14mu}{ratio}\mspace{14mu}{r\left( {{\overset{\sim}{\beta}}_{g}^{(i)},\beta_{g}^{({i - 1})}} \right)}} = \frac{{f\left( {{{\overset{\sim}{\beta}}_{g}^{(i)}❘\gamma_{g}^{(i)}},\kappa^{(i)},x^{(i)},y} \right)}{q\left( {\beta_{g}^{({i - 1})}❘{\overset{\sim}{\beta}}_{g}^{(i)}} \right)}}{{f\left( {{\beta_{g}^{({i - 1})}❘\gamma_{g}^{(i)}},\kappa^{(i)},x^{(i)},y} \right)}{q\left( {{\overset{\sim}{\beta}}_{g}^{(i)}❘\beta_{g}^{({i - 1})}} \right)}}$and accept the proposed candidate with the probability min{1, r({tildeover (β)}_(g) ^((i)), β_(g) ^((i−1)))}. 8: end for 9: Simulate κ^((i))as follows: Generate {tilde over (κ)}^((i)) ~ q(κ|κ^((i−1))).${{Compute}\mspace{14mu}{the}\mspace{14mu}{ratio}\mspace{14mu}{r\left( {{\overset{\sim}{\kappa}}^{(i)},\kappa^{({i - 1})}} \right)}} = \frac{{f\left( {{{\overset{\sim}{\kappa}}^{(i)}❘\gamma_{1}^{(i)}},\ldots\mspace{14mu},\gamma_{G}^{(i)},\beta_{1}^{(i)},\ldots\mspace{14mu},\beta_{G}^{(i)},x^{(i)},y} \right)}{q\left( {\kappa^{({i - 1})}❘{\overset{\sim}{\kappa}}^{(i)}} \right)}}{{f\left( {{\kappa^{({i - 1})}❘\gamma_{1}^{(i)}},\ldots\mspace{14mu},\gamma_{G}^{(i)},\beta_{1}^{(i)},\ldots\mspace{14mu},\beta_{G}^{(i)},x^{(i)},y} \right)}{q\left( {{\overset{\sim}{\kappa}}^{(i)}❘\kappa^{({i - 1})}} \right)}}$and accept the proposed candidate with the probability min{1, r({tildeover (κ)}^((i)), κ^((i−1)))}. 10: Set i ← i + 1. 11: Until Convergence

The technical results of different embodiments of the inventive methodswill now be discussed with reference to FIG. 2-20.

The anatomical results reported in FIGS. 2-10 and 13 were obtained inthe following conditions. Experiments have been conducted on real datasets comprising 256×256×14 Gradient-Echo (GE) anatomical with0.93×0.93×8 mm3 spatial resolution. GE anatomical images were acquiredwith TE/TR=10/500 ms and BW=31.25 kHz. Note also that these images havebeen acquired using acceleration factors R=2 and R=4 on a Signa 1.5Tesla GE Healthcare scanner with an eight-channel head coil.Interestingly, the scanning time of anatomical data lasted 5˜mn innon-parallel imaging, while acquisition duration was decreased to 3 mn10 s and 2 mn 20 s in parallel imaging with R=2 and R=4, respectively.

FIGS. 2 and 3 show 9 anatomical slices of a human brain obtained usingSENSE with Tikhonov regularization. The Tikhonov parameter has beenfixed manually. The reduction factor is R=2 for FIG. 2, and R=4 for FIG.3. Despite the regularization, some aliasing artifacts are visible inthe form of curved lines, particularly on FIG. 3. Some oversmoothing ofthe images is also visible.

FIG. 4 has been obtained, in the same experimental conditions, using TVregularization (left panel: R=2; right panel: R=4—a single slice isshown) with a manually tuned regularization parameter κ. Quite strongaliasing artifacts and some “staircase” defects are visible; thesedefects can be mitigated by increasing the value of κ, at the expense ofthe informational content of the images.

FIGS. 5 and 6 have been obtained, in the same experimental conditions,using an autocalibrated 2D wavelet transform-based regularizationscheme, as described above; this method will be called 2D-UWR-SENSE,where “UWR” stands for “Unconstrained Wavelet Regularization”. Again,R=2 for FIG. 5 and R=4 for FIG. 6.

More precisely, dyadic (M=2) Symmlet orthonormal wavelet basesassociated with filters of length 8 were used over j_(max)=3 resolutionlevels. Regarding the wavelet coefficients, the prior of equation (10)has been employed.

The related hyper-parameters (a couple of hyper-parameters is fitted forreal/imaginary parts of each subband, i.e. each approximation/detailcoefficients at each resolution level and orientation) were estimatedusing the Bayesian approach described above. Full FOV imagereconstruction was then performed using wavelet-based regularization.

The smoothing effects observed in FIGS. 2 and 3 with Tikhonovregularization no longer exist in the WT regularized images of FIGS. 5and 6, where a quite accurate reconstruction is performed within thebrain mask without introducing the staircase effects observed with TVregularization (see FIG. 4).

A further improvement of the image quality can be obtained byincorporating an additional constraint in the method described hereabovein order to better regularize artifact regions. The ensuing algorithm iscalled CWR-SENSE, where “CWR” stands for “Constrained WaveletRegularization”.

This method implies imposing local lower and upper bounds on the imageintensity values in artifact areas, regardless of their shape and/orlocation. These bounds define the nonempty closed convex set:C={ρ∈

^(K) |∀r∈{1, . . . ,Y/R}×{1, . . . ,X},ρ(r)∈C _(r)}  (53)

where the constraint introduced on the range values at position r∈2{1, .. . , Y/R}×{1, . . . , X} is modeled by:C _(r) ={ξ∈

|Re(ξ)∈

_(r) ^(Re) ,Im(ξ)∈

_(r) ^(Im)}  (54)with

_(r) ^(Re)=[I_(min,r) ^(Re), I_(max,r) ^(Re)] and

_(r) ^(Im)=[I_(min,r) ^(Im), I_(max,r) ^(Im)]. The optimality criterionbecomes then:

_(CWT)(ζ)=

_(WT)(ζ)+i _(C)·(ζ)  (55)where C*={ζ∈

^(K)|T*ζ∈C} and is the indicator function of the closed convex seti_(C)* defined by:

$\begin{matrix}{{\forall{\zeta \in {\mathbb{C}}^{K}}},\mspace{14mu}{{i_{C^{*}}(\zeta)} = \left\{ \begin{matrix}0 & {{{if}\mspace{14mu}\zeta} \in C^{*}} \\{+ \infty} & {{otherwise}.}\end{matrix} \right.}} & (56)\end{matrix}$

A morphological gradient [Serra, 1982] can used to detect artifactregions (very low/high transitions in the gradient image) on which theadditional convex constraints are applied. The upper and lower boundsthat define the convex sets C_(r) are spatially varying since theydepend on r. They can thus be computed using a morphological opening andclosing operations applied to the basic-SENSE reconstructed image(auxiliary, or reference, image) in order to discard very low and highintensities.

FIGS. 7 and 8 illustrate the results of the CWR-SENSE method, for R=2and R=4 respectively. It can be seen that the surviving artifacts inFIGS. 5 and 6 have now been removed due to the anisotropic smoothingusing the additional convex constraint.

From a quantitative point of view, significant improvements wereachieved by the UWR/CWR-SENSE algorithms in comparison with basic-SENSEand Tikhonov reconstructions. The following table 1 reports thesignal-noise-ratio (SNR) values in dB corresponding to the basic-SENSE,Tikhonov regularization and the proposed UWR/CWR-SENSE techniques forthe illustrated slices of the anatomical brain volume shown in FIGS. 3,6 and 8 (R=4). On average, the gain drew from the proposed constrainedregularization strategy amounts to 1.08 dB and a better visual quality.

TABLE 1 SENSE Tikhonov UWR-SENSE CWR-SENSE Slice #1 14.36 14.63 14.7714.95 Slice #2 11.55 11.62 12.01 12.53 Slice #3 12.95 13.44 14.02 14.22Slice #4 9.24 9.48 10.01 10.30 Slice #5 11.50 11.79 12.06 12.25 Slice #69.68 9.87 10.12 10.32 Slice #7 11.00 11.56 11.85 12.00 Slice #8 12.1612.49 13.00 13.36 Slice #9 13.78 14.77 15.83 16.04 Volume average 11.8012.18 12.63 12.88

The influence of the choice of the wavelet basis has also been studied.More particularly, four different bases have been considered: dyadicSymmlet 8, dyadic Daubechies 8, dyadic Haar and Meyer with M=4 bands[Chaux et al., 2006b]. The first three bases give quite similar results(apart from some blocking effects which only occur with the Haar basis),the dyadic Symmlet 8 leading to a slightly higher SNR. On the contrary,the Meyer 4-band wavelet basis leads to a significantly lower SNR.

FIG. 9 shows the results obtained by combined wavelet-total variationregularization (algorithm 2) using a dyadic (M=2) Symmlet orthonormalwavelet basis associated with filters of length 8 over j_(max)=3resolution levels, with a reduction factor R=4. FIG. 10 show the resultsobtained by combined wavelet-total variation regularization using aredundant wavelet frame constituted by a Union of 2 Orthonormal Bases(U2OB) with Symmlet 8 and Symmlet 4 filters. The hyper-parameters of thewavelet prior and the TV regularization parameter were estimated usingthe approach described above. These figures show that the reconstructedimages present better regularity than the ones reconstructed using theUWR-SENSE algorithm (FIG. 6). However, from a visual viewpoint, usingredundant WT do not necessarily lead to better reconstruction quality.From a quantitative viewpoint, SNR values in dB are provided in thetable 2 below. Comparisons with SNR values in table 1 confirm theusefulness of combining WT and TV in a joint regularization framework.This table shows that a slight improvement of 0.02 dB is obtainedcompared to the UWR-SENSE algorithm when using Symmlet wavelets.However, more significant improvement (0.24 dB) is reached when usingthe U2OB redundant wavelet frame. It turns out then that even if thereconstruction performance obtained using an orthonormal wavelet basisand a redundant wavelet frame seems equivalent, SNR values indicate thatusing redundant frames is fruitful.

FIG. 13 shows the reconstructed images using TV regularization (left),3D-UWR-SENSE with a wavelet frame constituted by a union of twoorthonormal wavelet basis (middle) and a hybrid frame-TV regularizationmethod (right). The top row shows whole images, the bottom row amagnified detail. The two orthonormal bases whose union constitutes thewavelet frame used for reconstruction (middle and right columns) are aDaubechies basis of length 4 and a shifted Daubechies basis of length 8.

Three resolution levels have been used, which means that G=20 groups ofwavelet coefficients are considered. For the frame representations, thehyperparameters have been chosen using the hybrid MCMC algorithmdescribed above.

TABLE 2 OB U2OB Slice #1 14.85 15.16 Slice #2 12.18 12.35 Slice #3 13.8514.32 Slice #4 9.78 9.85 Slice #5 12.17 12.62 Slice #6 10.20 10.17 Slice#7 11.75 12.04 Slice #8 13.15 13.44 Slice #9 15.94 15.88 Volume average12.65 12.87

The anatomical results reported in FIGS. 11 and 12. were obtained in thefollowing conditions. The anatomical MRI scan was performed on a 3T TimTrio Siemens scanner using a 3D T1-weighted MP-RAGE pulse sequence and amatrix array head coil consisting of 32 receive channels. The scanparameter were chosen as follows: slice orientation=sagittal(Right->Left), slice thickness=1.1 mm; TE=2.98 ms, TR=2300 ms; TI=900ms, Flip Angle: 9°, BW=61 kHz, single shot. The field of View was256×240×176 mm3 with a matrix size of 256×240×160 corresponding toanisotropic resolution of 1×1×1.1 mm3. The scanning time was TA_conv=9min 14 s using conventional MRI i.e. without acceleration (neither pMRInor partial Fourier). The use of a 6/8 partial Fourier scheme enabled todecrease the scanning time to 7 min 46 s. The parallel MRI data werecollected without partial Fourier using an acceleration factor R=2 orR=4. The respective scanning times for these exams were 5 min 03 s and 2min 59 s and not TA_conv/R as expected. The reason for this lies in theSiemens k-space sampling strategy. The manufacturer adopts a fullsampling scheme for the 24 central lines inducing an actual accelerationfactor lower than the prescribed one (1.83 instead of 2 and 3.09 insteadof 4).

FIG. 11-12 illustrate different MRI reconstruction algorithms fromsingle subject T1-weighted MRI data acquired at 3 Tesla (Siemens TimTrio) using a matrix array 32-channel head coil. FIGS. 11 and 12 showaxial and coronal views, respectively for the same subject. Top row inthese figures shows the ground truth i.e. the reconstruction performedfrom full k-space acquisition. Then, the mSENSE (middle row) and3D-UWR-SENSE (bottom row) algorithms are compared in a noisy context iefor R=4. For the latter algorithm the results have been obtained for J=2resolution levels and Daubechies wavelets (non-redundant). Also, thehyperparameters of the wavelet representation have been estimated usinga complete-data maximum likelihood procedure, as detailed above in themanuscript. MSENSE reconstruction artefacts appearing as white matterspots on the grey matter's boundaries are indicated by white circles inboth views outlines. Clearly, the 3D-UWR-SENSE algorithm outperforms themSENSE technique since these artefacts do not appear at the samelocation. Also, the elliptical-shaped artefacts in the centre of theaxial view are strongly filtered by the inventive wavelet-basedalgorithm. Coronal views confirm the substantial noise reduction in thebrain stem and subcortical regions using the proposed invention.

The results obtained using 2D wavelet bases are not shown here, as theyare qualitatively close to the ones reported in FIGS. 11 and 12.However, the 3D wavelet-based reconstruction quantitatively outperformsthe 2D-wavelet-based reconstruction in terms of Signal to Noise Ratio(SNR). Note that the 3D reconstruction generated a SNR improvement of1.3 dB with respect to 2D technique.

Until here, only results regarding the case of static, anatomicalimaging have been discussed. However, as discussed above, the presentinvention also applies to the space-time regularization of 4D fMRI imageseries, obtained e.g. by echoplanar imaging (EPI). The correspondingresults are illustrated by FIGS. 14-20.

For validation purpose, fMRI data were acquired on a 3 T Siemens Triomagnet using a Gradient-Echo EPI (GE-EPI) sequence (TE=30 ms, TR=2400ms, slice thickness=3 mm, transversal orientation, FOV=192 mm²) during acognitive localizer [Pinel et al, 2007] protocol. This experiment hasbeen designed to map auditory, visual and motor brain functions as wellas higher cognitive tasks such as number processing and languagecomprehension (listening and reading). It consisted of a single sessionof N_(r)=128 scans. The paradigm was a fast event-related designcomprising sixty auditory, visual and motor stimuli, defined in tenexperimental conditions (auditory and visual sentences, auditory andvisual calculations, left/right auditory and visual clicks, horizontaland vertical checkerboards). An L=32 channel coil was used to enableparallel imaging.

For each of 15 subjects, fMRI data were collected at the 2×2 mm² spatialin-plane resolution using different reduction factors (R=2 or R=4).Based on the raw data files delivered by the scanner, reduced FOV EPIimages were reconstructed using two specific treatments:

i) k-space regridding to account for the non-uniform k-space samplingduring readout gradient ramp, which occurs in fast MRI sequences likeGE-EPI;

ii) Nyquist ghosting correction to remove the odd-even echoinconsistencies during k-space acquisition of EPI images.

Then, the 4D-UWR-SENSE and the already discussed 3D-UWR-SENSE (simplycalled “UWR-SENSE” for short) algorithms have been utilized in a finalstep to reconstruct the full FOV EPI images and compared to the mSENSEsolution (mSENSE is the unregularized SENSE algorithm implemented bySiemens scanners).

FIG. 14 compares the two pMRI reconstruction algorithms to illustrate onaxial, coronal and sagittal slices how the mSENSE reconstructionartifacts have been removed using the 4D-UWR-SENSE approach. The mSENSEreconstructed images actually present large artifacts located both atthe center and boundaries of the brain in sensory and cognitive regions(temporal lobes, frontal and motor cortices); on the figure, theartifacts are outlined by ellipses. This results in SNR loss and thusmay have a dramatic impact for activation detection in these brainregions.

Irrespective of the reconstruction pipeline, the full FOV fMRI imageswere then preprocessed using the SPM5 software(http://www.fil.ion.ucl.ac.uk/spm/software/spm5/): preprocessinginvolved realignment, correction for motion and differences in sliceacquisition time, spatial normalization, and smoothing with an isotropicGaussian kernel of 4 mm full-width at half-maximum. Anatomicalnormalization to MNI space was performed by coregistration of thefunctional images with the anatomical T1 scan acquired with the thirtytwo channel-head coil. Parameters for the normalization to MNI spacewere estimated by normalizing this scan to the T1 MNI template providedby SPM5, and were subsequently applied to all functional images.

For performing subject-level analysis, a General Linear Model (GLM) wasconstructed to capture stimulus-related BOLD response. The design matrixrelied on ten experimental conditions and thus made up of twenty-oneregressors corresponding to stick functions convolved with the canonicalHemodynamic Response Function (HRF) and its first temporal derivative,the last regressor modeling the baseline. This GLM was then fitted tothe same acquired images and reconstructed using the Siemensreconstructor, UWR-SENSE and 4D-UWR-SENSE.

Contrast estimate images for motor responses and higher cognitivefunctions (computation, language) were subjected to further analyses atthe subject and group levels. These two contrasts are complementarysince the expected activations lie in different brain regions and thuscan be differentially corrupted by reconstruction artifacts.

More precisely, were studied:

-   -   the Auditory computation vs. Auditory sentence (aC-aS) contrast        which is supposed to elicit evoked activity in the frontal and        parietal lobes, since solving mental arithmetic task involves        working memory and more specifically the intra-parietal sulcus;    -   the Left click vs. Right click (Lc-Rc) contrast for which we        expect evoked activity in the right motor cortex (precentral        gyrus, middle frontal gyrus). Indeed, the Lc-Rc contrast defines        a compound comparison which involves two motor stimuli which are        presented either in the visual or auditory modality. This        comparison aims therefore at detecting lateralization effect in        the motor cortex.

These two contrasts were chosen because they summarized well differentsituations (large vs small activation clusters, distributed vs focalactivation pattern, bilateral vs unilateral activity) that occurred forthis paradigm when looking at sensory areas (visual, auditory, motor) orregions involved in higher cognitive functions (reading, calculation).

FIG. 15 shows subject-level student-t maps superimposed to anatomicalMRI for the aC-aS contrast. Data have been reconstructed using themSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively, with R=2 (top of thefigure) and R=4 (bottom of the figure). The neurological convention(left is left) is adopted. The cross shows the maximum activation peak.

For the most significant slice and R=2, all pMRI reconstructionalgorithms succeed in finding evoked activity in the left parietal andfrontal cortices, more precisely in the inferior parietal lobule andmiddle frontal gyrus. However, for R=4 only UWR-SENSE and4D-UWR-SENSE—and preferentially the latter—enable to retrieve reliablefrontal activity elicited by mental calculation, which is lost by themSENSE algorithm. From a quantitative viewpoint, the proposed4D-UWR-SENSE algorithm finds larger clusters whose local maxima are moresignificant than the ones obtained using mSENSE and UWR-SENSE, asreported in Table 3. Concerning the most significant cluster for R=2,the peak positions remain stable whatever the reconstruction algorithm.However, examining their significance level, one can first measure thebenefits of wavelet-based regularization when comparing UWR-SENSE withmSENSE results and then additional positive effects of temporalregularization and 3D wavelet decomposition when looking at the4D-UWR-SENSE results. These benefits are also demonstrated for R=4.Table 3 shows the significant statistical results at the subject-levelfor the aC-aS contrast (corrected for multiple comparisons at thesignificance level of α=0.05, which means that the null hypothesis isrejected if the p-value is smaller or equal to α).

The ‘p-value” of the statistical significance test should not beconfused with the “shape parameter” used in temporal regularization, seeequations 21 and 22, or with the parameter “p” used to illustrate theMCMC algorithm.

TABLE 3 cluster-level voxel-level p-value Size p-value T-score PositionR = 2 mSENSE <10⁻³ 320 <10⁻³ 6.40 −32 −76 45 <10⁻³ 163 <10⁻³ 5.96 −4 −7054 <10⁻³ 121 <10⁻³ 6.34 34 −74 39 <10⁻³ 94 <10⁻³ 6.83 −38 4 24 UWR-SENSE<10⁻³ 407 <10⁻³ 6.59 −32 −76 45 <10⁻³ 164 <10⁻³ 5.69 −6 −70 54 <10⁻³ 159<10⁻³ 5.84 32 −70 39 <10⁻³ 155 <10⁻³ 6.87 −44 4 24 4D-UWR- <10⁻³ 454<10⁻³ 6.54 −32 −76 45 SENSE <10⁻³ 199 <10⁻³ 5.43 −6 26 21 <10⁻³ 183<10⁻³ 5.89 32 −70 39 <10⁻³ 170 <10⁻³ 6.90 −44 4 24 R = 4 mSENSE <10⁻³ 580.028 5.16 −30 −72 48 4D-UWR- <10⁻³ 94 0.003 5.91 −32 −70 48 SENSE <10⁻³60 0.044 4.42 −6 −72 54 4D-UWR- <10⁻³ 152 <10⁻³ 6.36 −32 −70 48 SENSE<10⁻³ 36 0.009 5.01 −4 −78 48 <10⁻³ 29 0.004 5.30 −34 6 27

FIG. 16 illustrates between-subject variability of detected activationfor the aC-aS contrast at R=2. Indeed, when comparing subject-levelstudent-t maps reconstructed using the different pipelines (R=2), it canbe observed that the mSENSE algorithm fails to detect any activationcluster in the expected regions for the second subject. In contrast, the4D-UWR-SENSE method retrieves more coherent activity while not exactlyat the same position as for the first subject.

FIG. 17 shows subject-level student-t maps superimposed to anatomicalMRI for the Lc-Rc contrast. Data have been reconstructed using themSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively.

It can be seen that all reconstruction methods enable to retrieveexpected activation in the right precentral gyrus. However, when lookingmore carefully at the statistical results (see Table 4), the UWR-SENSEand more preferentially the 4D-UWR-SENSE algorithms retrieve anadditional cluster in the right middle frontal gyrus. On data acquiredwith R=4, the same Lc-Rc contrast elicits similar activations, i.e. inthe same region. As it can be seen on the bottom of the figure, thisactivity is enhanced when pMRI reconstruction is performed accordingwith the methods of the invention.

Quantitative results in Table 4 confirm numerically what can be observedin the figure: larger clusters with higher local t-scores are detectedusing the 4D-UWR-SENSE algorithm, both for R=2 and R=4. More precisely,table 4 shows significant statistical results at the subject-level forthe Lc-Rc contrast (corrected for multiple comparisons at thesignificance level of α=0.05, which means that the null hypothesis isrejected if the p-value is smaller or equal to α).

TABLE 4 cluster-level voxel-level p-value Size p-value T-score PositionR = 2 mSENSE <10⁻³ 79 <10⁻³ 6.49 38 −26 66 UWR-SENSE <10⁻³ 144 0.0045.82 40 −22 63 0.03  21 0.064 4.19 24 −8 63 4D-UWR- <10⁻³ 172 0.001 6.7834 −24 69 SENSE <10⁻³ 79 0.001 6.49 38 −26 66 R = 4 mSENSE 0.006 210.295 4.82 34 −28 63 UWR-SENSE <10⁻³ 33 0.120 5.06 40 −24 66 4D-UWR-<10⁻³ 51 0.006 5.57 40 −24 66 SENSE

FIG. 18 reports on the robustness of the proposed pMRI pipeline to thebetween-subject variability for this motor contrast. Since sensoryfunctions are expected to generate larger BOLD effects (higher SNR) andappears more stable, our comparison takes place at R=4. Twosubject-level student-t maps reconstructed using the different pMRIalgorithms are compared. For the second subject, one can observe thatthe mSENSE algorithm fails to detect any activation cluster in the rightmotor cortex. In contrast, the 4D-UWR-SENSE method retrieves morecoherent activity for this second subject in the expected region.

To summarize, on these two contrasts the 4D-UWR-SENSE algorithm alwaysoutperforms the alternative reconstruction methods in terms ofstatistical significance (number of clusters, cluster extent, peakvalues, . . . ) but also in terms of robustness.

Due to between-subject anatomical and functional variability,group-level analysis is necessary in order to derive robust andreproducible conclusions at the population level. For this validation,random effect analyses (RFX) involving fifteen healthy subjects havebeen conducted on the contrast maps we previously investigated at thesubject level. More precisely, one-sample Student-t test was performedon the subject-level contrast images (e.g., Lc-Rc, aC-aS, images) usingSPM5.

FIG. 19 shows group-level student-t maps for the aC-aS contrast wheredata have been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSEfor R=2 and R=4. Neurological convention has been used. Arrows indicatethe global maximum activation peak.

These maps illustrate that irrespective of the reconstruction methodlarger and more significant activations are found on datasets acquiredwith R=2 given the better SNR. Second, for R=2, visual inspectionconfirms that only the 4D-UWR-SENSE algorithm allows to retrievesignificant bilateral activations in the parietal cortices (see axialMIP slices) in addition to larger cluster extent and a gain insignificance level for the stable clusters across the differentreconstructors. Similar conclusions can be drawn when looking at thebottom of the figure, for R=4. Complementary results are available inTable 5 for R=2 and R=4 and numerically confirms this visual comparison:

-   -   Whatever the reconstruction method in use, the statistical        performance is much more significant using R=2, especially at        the cluster level since the cluster extent decreases by one        order of magnitude.    -   Voxel and cluster-level results are enhanced using the        4D-UWR-SENSE approach instead of the mSENSE or UWR-SENSE.

TABLE 5 Significant statistical results at the group-level for the aC-aScontrast (corrected for multiple comparisons at p = 0.05). cluster-levelvoxel-level p-value Size p-value T-score Position R - 2 mSENSE <10⁻³ 3610.014 7.68 −6 −22 45 <10⁻³ 331 0.014 8.23 −40 −38 42 <10⁻³ 70 0.014 7.84−44 6 27 UWR-SENSE <10⁻³ 361 0.014 7.68 −6 22 45 <10⁻³ 331 0.014 7.68−44 −38 42 <10⁻³ 70 0.014 7.84 −44 6 27 4D-UWR- <10⁻³ 441 <10⁻³ 9.45 −32−50 45 SENSE <10⁻³ 338 <10⁻³ 9.37 −6 12 45 <10⁻³ 152 0.010 7.19 30 −6448 R - 4 mSENSE 0.003 14 0.737 5.13 −38 −42 51 UWR-SENSE <10⁻³ 41 0.2745.78 −50 −38 −48 <10⁻³ 32 0.274 5.91 2 12 54 4D-UWR- <10⁻³ 37 0.268 6.46−40 −40 54 SENSE <10⁻³ 25 0.268 6.37 −38 −42 36 <10⁻³ 18 0.273 5 −42 836

FIG. 20 reports similar group-level MIP results for R=2 and R=4concerning the Lc-Rc contrast. It is shown that whatever theacceleration factor R in use, our pipeline enables to detect a much morespatially extended activation area in the motor cortex. This visualinspection is quantitatively confirmed in Table 6 when comparing thedetected clusters using the 4D-UWR-SENSE approach with those found bymSENSE, again irrespective of R. Finally, the 4D-UWR-SENSE algorithmoutperforms the UWR-SENSE one, which corroborates the benefits of theproposed spatio-temporal regularization scheme.

TABLE 6 Significant statistical results at the group-level for the Lc-Rccontrast (corrected for multiple comparisons at p = 0.05). cluster-levelvoxel-level p-value Size p-value T-score Position R = 2 mSENSE <10⁻³ 354<10⁻³ 9.48 38 −22 54 0.001 44 0.665 6.09 −4 −68 −24 UWR-SENSE <10⁻³ 3500.005 9.83 36 −22 57 <10⁻³ 35 0.286 7.02 4 −12 51 4D-UWR- <10⁻³ 3770.001 11.34 36 −22 57 SENSE <10⁻³ 53 <10⁻³ 7.50 8 −14 51 <10⁻³ 47 <10⁻³7.24 −18 −54 −18 R = 4 mSENSE <10⁻³ 38 0.990 5.97 32 −20 45 UWR-SENSE<10⁻³ 163 0.128 7.51 46 −18 60 4D-UWR- <10⁻³ 180 0.111 7.61 46 −18 60SENSE

REFERENCES

-   [Bertsekas, 1995] Bertsekas, D. P. (1995). Nonlinear programming,    Second Edition. Athena Scientific, Belmont, USA. In particular,    pages 159-165.-   [Block et al., 2007] Block, K. T., Uecker, M., and Frahm, J. (2007).    Undersampled radial MRI with multiple coils. Iterative image    reconstruction using a total variation constraint. Magnetic    Resonance in Medicine, 56(7):1086-1098.-   [Chaâri et aI.2008]: L. Chaâri, J.-C. Pesquet, A.    Benazza-Benyahia, P. Ciuciu, Autocalibrated Parallel MRI    Reconstruction in the Wavelet Domain, in: IEEE International    Symposium on Biomedical Imaging, Paris, France, 2008, pp. 756-759.-   [Chaâri et al. 2009]: L. Chaâri, J.-C. Pesquet, Ph. Ciuciu, and A.    Benazza-Benyahia An Iterative Method for Parallel MRI SENSE-based    Reconstruction in the Wavelet Domain., arXiv:0909.0368v1 [math.OC]-   [Chaux et al., 2007]: Chaux, C., Combettes, P., Pesquet, J.-C., and    Wajs, V. (2007). A variational formulation for frame-based inverse    problems. Inverse Problems, 23(4):1495-1518.-   [Combettes and Pesquet, 2008]: Combettes, P. L. and Pesquet, J.-C.    (2008). A proximal decomposition method for solving convex    variational inverse problems. Inverse Problems, 24(6):27.-   [Daubechies et al., 2004]: I. Daubechies, M. Defrise, C. DeMol, An    iterative thresholding algorithm for linear inverse problems with a    sparsity constraint, Communications on Pure and Applied Mathematics    57 (11) (2004) 1413-1457.-   [Dempster, 1997]: Dempster, A., Laird, A., and Rubin, D. (1977).    Maximum likelihood from incomplete data via the EM algorithm (with    discussion). Journal of the Royal Statistical Society, Series B,    39:1-38.-   [Dobigeon et al, 2007]: Dobigeon, N. and Tourneret, J.-Y. (2007).    Truncated multivariate Gaussian distribution on a simplex. Technical    report, University of Toulouse.-   [Donoho, 1995]: D. Donoho. De-noising by soft-thresholding. IEEE    Transactions on Information Theory, 41(3):613-627, 1995.-   [Geman, 1984]: S. Geman and D. Geman, “Stochastic relaxation, Gibbs    distribution and the Bayesian restoration of image,” IEEE Trans.    Pattern Anal. Mach. Intell., vol. 6, pp. 721-741, 1984.-   [Griswold et al., 2002]: M. A. Griswold, P. M. Jakob, R. M.    Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, A. Haase,    Generalized autocalibrating partially parallel acquisitions GRAPPA,    Magnetic Resonance in Medicine 47 (6) (2002) 1202-1210.-   [Gupta et al., 1997]: Gupta, A. K. and Song, D. (1997). Journal of    Statistical Planning and Inference. Volume 60, Pages 241-260.-   [Hoge et al., 2005]: W. S. Hoge, D. H. Brooks, B. Madore, W. E.    Kyriakos, A tour of accelerated parallel MR imaging from a linear    systems perspective, Concepts in Magnetic Resonance 27A (1) (2005)    17-37.-   [Pesquet-Popescu, Pesquet] Ondelettes et applications, Techniques de    l'Ingénieur, traité Télécoms, TE 5 215-   [Pinel et al., 2007]: P. Pinel, B. Thirion, S. Meriaux, A.    Jobert, J. Serres, D. Le Bihan, J.-B. Poline, and S. Dehaene, “Fast    reproducible identification and large-scale databasing of individual    functional cognitive networks” BMC Neurosci., vol. 8, no. 1, pp. 91,    October 2007.-   [Pruessmann et al., 1999]: K. P. Pruessmann, M. Weiger, M. B.    Scheidegger, P. Boesiger, SENSE: sensitivity encoding for fast MRI,    Magnetic Resonance in Medicine 42 (5) (1999) 952-962.-   [Raj et al., 2007]: Raj, A., Singh, G., Zabih, R., Kressler, B.,    Wang, Y., Schuff, N., and Weiner, M. (2007). Bayesian parallel    imaging with edge-preserving priors. Magnetic Resonance in Medicine,    57(1):8-21.-   [Robert, 1995]: “Simulation of truncated normal variables,”    Statistics and Computing, vol. 5, pp. 121-125,1995.-   [Serra, 1982] Serra, J. (1982). Image Analysis and Mathematical    Morphology. Academic Press, London.-   [Sodickson et al., 1997]: Sodickson, D. K. and Manning, W. J.    (1997). Simultaneous acquisition of spatial harmonics (SMASH): fast    imaging with radiofrequency coil arrays. Magnetic Resonance in    Medicine, 38(4):591-603.-   [Ying et al., 2004]: L. Ying, D. Xu, Z.-P. Liang, On Tikhonov    Regularization for image reconstruction in parallel MRI, in: IEEE    Engineering in Medicine and Biology Society, San Francisco, USA,    2004, pp. 1056-1059.

The invention claimed is:
 1. A method of performing parallel magneticresonance imaging of a body in a magnetic resonance imaging system,comprising: acquiring a set of elementary magnetic resonance images ofsaid body from respective radio frequency (RF) receiving antennas havingknown or estimated sensitivity maps and noise covariance matrices, saidacquired elementary magnetic resonance images being under-sampled ink-space; and performing with a processor a regularized reconstruction oneach acquired magnetic resonance image of said body; where said step ofperforming regularized reconstruction on a magnetic resonance image ofsaid body is carried out in a discrete frame space by minimizing a costfunction, with the cost function comprising: an error term, thatrepresents a statistical likelihood of a reconstructed image, obtainedfrom said acquired elementary magnetic resonance images; and a framepenalty term, representative of a deviation between an actualstatistical distribution of frame coefficients of said reconstructedimage of said body and an a priori distribution of said framecoefficients; said a priori distributions of the frame coefficients ofthe reconstructed image of said body being estimated on the basis of anauxiliary magnetic resonance image of said body; and wherein saidauxiliary magnetic resonance image of said body is reconstructed by saidprocessor from said acquired elementary magnetic resonance images, usingan algorithm chosen between: a constrained autocalibrated 2D wavelettransform-based regularization scheme; an unconstrained autocalibrated2D wavelet transform-based regularization scheme; an unconstrainedautocalibrated 3D wavelet transform-based regularization scheme; and anunconstrained autocalibrated 4D wavelet transform-based regularizationscheme.
 2. A method according to claim 1, wherein said error term isrepresentative of a statistical neg-log-likelihood of said reconstructedimage of said body.
 3. A method according to claim 1, wherein said stepof performing regularized reconstruction of a magnetic resonance imageof said body is carried out by maximizing, in said frame space, a fullposterior distribution of a set of frame coefficients that are utilizedin defining an image of the body, from said acquired elementary magneticresonance images and said a priori distribution of the framecoefficients.
 4. A method according to claim 1, wherein said auxiliarymagnetic resonance image of said body is reconstructed from saidacquired elementary magnetic resonance images by using a SENSitivityEncoding-SENSE-algorithm.
 5. A method according to claim 4 wherein saidauxiliary magnetic resonance image of said body is reconstructed fromsaid acquired elementary magnetic resonance images by using an algorithmchosen between: an unregularized SENSE algorithm; a SENSE algorithmregularized in image space; and a SENSE algorithm regularized ink-space.
 6. A method according to claim 1 wherein a generalizedGauss-Laplace a priori statistical distribution of said framecoefficients is assumed, and parameters of said Gauss-Laplace a prioristatistical distribution are estimated on the basis of said auxiliarymagnetic resonance image of said body, by using either a statisticalmaximum-likelihood or a posterior mean estimator.
 7. A method accordingto claim 1, wherein said error term that represents a statisticallikelihood of a reconstructed image, is a quadratic mean error term. 8.A method according to claim 1, wherein said acquired elementary magneticresonance images of the body are three-dimensional images, and whereinsaid step of performing with a processor said regularized reconstructionof a magnetic resonance image is carried out in a discretethree-dimensional frame space.
 9. A method according to claim 8, whereinsaid acquired three-dimensional elementary magnetic resonance images ofthe body are obtained by stacking bi-dimensional elementary magneticresonance images of slices of the object to be imaged.
 10. A methodaccording to claim 1 wherein said step of performing with a processor,said regularized reconstruction of a magnetic resonance image is basedon a redundant wavelet frame representation.
 11. A method according toclaim 1 wherein said step of performing regularized reconstruction of amagnetic resonance is based on a non-redundant wavelet regularizationscheme.
 12. A method according to claim 1 wherein said cost functionalso comprises at least one spatial domain penalty term chosen between:a total variation norm of the reconstructed image; and a convexconstraint of the reconstructed image.
 13. A method of performingdynamic parallel magnetic resonance imaging of a body, comprising:acquiring a set of time series of elementary magnetic resonance imagesof said body from respective radio frequency (RF) receiving antennashaving known or estimated sensitivity maps and noise covariancematrices, said acquired elementary images being under-sampled ink-space; and performing regularized reconstruction of a time series ofmagnetic resonance images of said body; where said step of performingwith a processor a regularized reconstruction on a time series ofelementary magnetic resonance images of said body being carried out byminimizing a non-differentiable cost function, with saidnon-differentiable cost function comprising: an error term, thatrepresents a statistical likelihood of each reconstructed image, giventhe acquired time series of corresponding acquired elementary magneticresonance images of said body; an edge-preserving temporal penalty term,that is representative of a pixel-by-pixel or voxel-by-voxel differencebetween consecutive regularized reconstructed images occurring in theseries, and wherein said time series of magnetic resonance image of saidbody is reconstructed by said processor from said acquired set of timeseries of elementary magnetic resonance images, using an algorithmchosen between: a constrained autocalibrated 2D wavelet transform-basedregularization scheme; an unconstrained autocalibrated 2D wavelettransform-based regularization scheme; an unconstrained autocalibrated3D wavelet transform-based regularization scheme; and an unconstrainedautocalibrated 4D wavelet transform-based regularization scheme.
 14. Amethod according to claim 13, wherein said temporal penalty term is aconvex edge-preserving function.
 15. A method according to claim 14,wherein said edge preserving temporal penalty term is a statistical Lpnorm with p≥1 and preferably 1≤p≤1.5.
 16. A method according to claim13, wherein said edge preserving temporal penalty term is given by thesum of a first partial temporal penalty term and a second partialtemporal penalty term in the non-differentiable cost function, wherein:the first partial temporal penalty term is representative ofpixel-by-pixel or voxel-by-voxel differences between each even-numberedimage of the series and a preceding odd-numbered image; and the secondpartial temporal penalty term is representative of pixel-by-pixel orvoxel-by-voxel differences between each odd-numbered image of the seriesand a preceding even-numbered image; with said non-differentiable costfunction being minimized by using proximity operators for said first andsecond partial temporal penalty terms that comprise the edge preservingtemporal penalty terms of the non-differentiable cost function.
 17. Amethod of performing dynamic parallel magnetic resonance imaging of abody according to claim 13, wherein said step of performing with saidprocessor a regularized reconstruction of a magnetic resonance image ofsaid body is carried out in a discrete frame space, and said costfunction also comprises a frame penalty term, representative of adeviation between statistical distributions of frame coefficientsobtained from each reconstructed image and said priori distributions ofsaid coefficients; said a priori distributions of the frame coefficientsof the reconstructed images being estimated from an auxiliary magneticresonance image of said body.
 18. A method according to claim 17,wherein said elementary magnetic resonance images of said body arethree-dimensional images, and wherein said discrete frame space is adiscrete three-dimensional frame space.
 19. A method according to claim13 further comprising a step of automatically determining a weightingparameter of said edge preserving temporal penalty term by using astatistical maximum-likelihood estimator.
 20. A method according toclaim 19 comprising estimating said weighting parameter of the edgepreserving temporal penalty term for each pixel or voxel, or for eachset of neighboring pixels or for each set of neighboring voxels, of theimage to be reconstructed.
 21. A method according to claim 19, whereinsaid edge preserving temporal penalty term is a statistical Lp norm andwherein said weighting parameter of the edge preserving temporal penaltyterm and the parameter p are jointly determined by using saidmaximum-likelihood estimator as the weighting parameter.
 22. A methodaccording to claim 21, wherein said weighting parameter of the edgepreserving temporal penalty term and the parameter p are jointlydetermined by using said statistical maximum-likelihood estimator underthe constraint p≥1.
 23. A method according to claim 13 wherein saiderror term depends on geometrical parameters defining a rigidtransformation of each of said acquired elementary magnetic resonanceimages of said body with respect to an additional elementary magneticresonance image of said body taken as a reference, and wherein said stepof performing with a processor a regularized reconstruction of a timeseries of elementary magnetic resonance images of said body is carriedout by minimizing said non-differentiable cost function also withrespect to said geometrical parameters.
 24. A method according to claim1, wherein said elementary images of said body are acquired by echoplanar imaging.
 25. A method according to claim 1, wherein saidelementary images of said body are under-sampled with a reduction factorhigher than or equal to
 4. 26. A method according to claim 13, whereinsaid elementary images of said body are acquired by echo planar imaging.27. A method according to claim 13, wherein said elementary images ofsaid body are under-sampled with a reduction factor higher than or equalto 4.